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ABSTRACT 

Relativistic jets are a natural outcome of some of the most violent and spectacular astrophys- 
ical phenomena, such as the core collapse of massive stars in gamma-ray bursts (GRBs) and 
the accretion onto supermassive black holes in active galactic nuclei (AGN). It is generally 
accepted that these jets are powered electromagnetically, by the magnetised rotation of a cen¬ 
tral compact object (a black hole or neutron star). However, how the jets produce the observed 
emission and survive the propagation for many orders of magnitude in distance without be¬ 
ing disrupted by current-driven non-axisymmetric instabilities is the subject of active debate. 
We carry out time-dependent 3D relativistic magnetohydrodynamic simulations of relativis¬ 
tic, Poynting flux dominated jets that propagate in a spherically-symmetric power-law density 
distribution. The jets are launched self-consistently by the rotation of a strongly magnetised 
central compact object. This determines the natural degree of azimuthal magnetic field wind¬ 
ing, a crucial factor that controls jet stability. We find that the jets are susceptible to two types 
of instability: (i) a global, external kink mode that grows on long time scales and causes the 
jets to bodily bend sideways. Whereas this mode does not cause jet disruption over the simu¬ 
lated distances, it substantially reduces jet propagation speed. We show, via an analytic model, 
that the growth of the external kink mode depends on the slope of the ambient medium den¬ 
sity profile. In flat density distributions characteristic of galactic cores, an AGN jet may stall, 
whereas in stellar envelopes the external kink weakens as the jet propagates outward; (ii) a lo¬ 
cal, internal kink mode that grows over short time scales and causes small-angle magnetic 
reconnection and conversion of about half of jet electromagnetic energy flux into heat. Based 
on the robustness and energetics of the internal kink mode, we suggest that this instability is 
the main dissipation mechanism responsible for powering GRB prompt emission. 
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1 INTRODUCTION 

Relativistic jets are ubiquitous among astrophysical systems that 
involve accretion onto compact objects, such as gamma-ray bursts 
(GRBs), active galactic nuclei (AGNs) and microquasars (see, e.g. 
a review by Levinson 2006). It is commonly agreed that these 
jets are powered electromagnetically, most likely by the winding 
of magnetic field lines that thread a rotating central compact ob¬ 
ject (Blandford & Znajek 1977; Komissarov 2001; Tchekhovskoy 
et al. 2010a). The winding generates a Poynting flux dominated 
outflow, which eventually becomes the jet, at the expense of the ro¬ 
tational energy of the central engine. Although significant progress 
has been made recently in understanding how jets are formed 
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magnetically (see, e.g., Tchekhovskoy 2015 for a review), their 
physics, most notably the stability properties and the energy dis¬ 
sipation mechanisms, is the subject of active debate. In the con¬ 
text of non-relativistic jets, the expectation is that magnetised jets 
are strongly unstable to current-driven instabilities. For instance, 
Moll et al. (2008); Moll (2009) show that non-relativistic jets read¬ 
ily develop current-driven, non-axisymmetric, kink (m = 1) modes 
whose properties depend on the conditions at the jet base. In the 
core collapse of a massive star, heavily mass-loaded jets, which 
are produced by the rotation of a protoneutron star (Burrows et al. 
2007), were found to be so strongly unstable to the kink instability 
that they were unable to penetrate the star (Mosta et al. 2014). 

If this instability were to extend also to the relativistic, highly 
magnetised jets this would have serious implications on the prop¬ 
erties of relativistic jets and the engines that launch them. In the 
context of long-duration GRBs (Woosley 1993; MacFadyen & 
Woosley 1999) it would imply that magnetised jets are unable to 
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break out of the star, a necessary condition to form a GRB. There¬ 
fore, it would mean that GRB jets have to be created unmagnetised , 
making it impossible to power them by the electromagnetic spin- 
down of the central object. However, the observed properties of 
relativistic jets suggest otherwise. For example, the high power ob¬ 
served in GRB and AGN jets significantly challenges the known 
non-magnetic energy extraction mechanisms available in these ob¬ 
jects (e.g. Phinney 1982; Kawanaka et al. 2013; Leng & Giannios 
2014). Moreover, a recently discovered correlation between the jet 
magnetic field strength and the accretion disc luminosity in AGN 
(Zamaninasab et al. 2014), and several surprising features that are 
seen in GRB lightcurves and can be naturally produced by mag¬ 
netic jets (Tchekhovskoy & Giannios 2014) strongly support the 
magnetic origin of the jet power. 

On the other, hand if the magnetic jets were mostly stable it 
would be difficult to explain the high energy emission radiated from 
them (see, e.g., McKinney & Blandford 2009; Narayan et al. 2009). 
In a stable jet about half of the energy remains locked in the mag¬ 
netic form, out to very large distances (Tchekhovskoy et al. 2010b). 
In such a case internal shocks, which are commonly invoked to ex¬ 
plain the observed high energy emission, are weak (Kennel & Coro- 
niti 1984) and cannot accelerate efficiently the radiating electrons 
(e.g., Mimica et al. 2010; Narayan et al. 2011). Various alternative 
dissipation mechanisms have been discussed in this context, includ¬ 
ing striped wind like magnetic field configurations, which are sus¬ 
ceptible to reconnection (e.g., Drenkhahn & Spruit 2002; Giannios 
& Spruit 2005; Metzger et al. 2011; McKinney & Uzdensky 2012), 
microphysical energy dissipation mechanisms (e.g., Beloborodov 
2010; Meszaros & Rees 2011), and effects from intermittent engine 
(e.g. Granot et al. 2011). To date, no universally-accepted mecha¬ 
nism that is capable of efficiently converting jet magnetic energy 
into radiation exists. The presence of a local current-driven insta¬ 
bility in a mildly unstable jet could be sufficient for triggering mag¬ 
netic dissipation and powering the observed GRB emission (Lyu- 
tikov & Blandford 2003; Giannios 2008, 2012). This motivates a 
focused, high-resolution study of jet stability in the context of core¬ 
collapse GRB jets. 

Current-driven instability is mainly a 3D effect and can be 
highly non-linear. Thus 3D simulations are needed to study it. The 
standard approach for simulating magnetic jets is sending loops of 
magnetic field into the computational domain at a fixed rate (e.g., 
Mignone et al. 2010; Guan et al. 2014). In this case, the strengths of 
the poloidal magnetic field component, B p (i.e., lying in the plane 
passing through the jet axis), and azimuthal component, B$ (per¬ 
pendicular to the plane), are arbitrary. However, jet stability de¬ 
pends sensitively on the ratio between the two (Appl et al. 2000; 
Narayan et al. 2009; Mizuno et al. 2012). Thus, the results of such 
jet “injection” simulations reflect a particular (arbitrary) choice 
for the injection boundary condition, making it difficult to inter¬ 
pret them. Lately, Bromberg et al. (2014) showed analytically that 
Poynting-dominated jets that form at the center of collapsing stars 
are at least marginally stable and can punch through the stellar en¬ 
velope without being disrupted by magnetic instability. However, 
in the absence of a full 3D numerical simulation, they had to as¬ 
sume a ratio between B$ and B p , and did not address the question 
of magnetic energy dissipation. 1 


1 They employed a linear analysis and assumed that the strengths of the 
toroidal and poloidal field components are comparable (in the fluid frame). 
Such an assumption is usually attributed to jets that propagate through a 
pre-evacuated funnel and is commonly used in the studies of jet stability 
(e.g., Narayan et al. 2009; McKinney & Blandford 2009). 
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Figure 1 . A schematic view of a headed jet (panel a) that propagates in an 
ambient medium, and a headless jet (panel b) that propagates in a preexist¬ 
ing evacuated funnel. As the headed jet runs into the ambient gas, it slows 
down, its toroidal magnetic flux accumulates and grows in strength, and 
its outermost parts form a head (shown in dark blue), a working surface at 
which the jet drills through the ambient gas and in which the magnetic field 
strength is enhanced. This causes the formation of a bow shock (shown with 
the dark green line), and the shocked ambient gas forms the cocoon (shown 
in green) that collimates the jet into a cigar-like shape. In contrast, a head¬ 
less jet has no ambient gas to push through and is free to propagate along a 
pre-evacuated funnel. It propagates faster than headed jets and assumes the 
shape dictated by the funnel and/or ambient pressure profile. Because head¬ 
less jets do not have to push through the ambient gas, their toroidal field is 
weaker and they are more stable to kink instabilities than headed jets. 


In nature, the strengths of the poloidal and toroidal magnetic 
field components are tightly connected. This is because the jets 
are produced by the rotation of magnetised compact objects, so 
B^ emerges from the winding of B p . The winding creates an out¬ 
ward Poynting flux. In the presence of an ambient medium, which 
confines the electromagnetic outflow, the magnetic tension of the 
azimuthal field builds up, until it focuses the outflow into twin po¬ 
lar collimated jets (Lynden-Bell 1996), as illustrated in Fig. 1(a). 
These jets start propagating once their pressure becomes high 
enough to push the ambient medium aside. As the jet propagates, it 
develops a slow-moving “head”: a working surface at which the jet 
drills through the ambient material and which is shown in Fig. 1(a) 
in dark blue. The head blocks the free expansion of the toroidal 
magnetic flux in the jet and keeps the jet toroidal magnetic pres¬ 
sure high. As a result, the jet pushes against the head with a greater 
force. We term this type of jets headed jets. The relative strength 
of toroidal and poloidal fields in the jet is therefore linked with 
the properties of the ambient medium that collimates it. Thus, any 
attempt to analyse jet stability should take into consideration the 
presence of the ambient medium and its effect on the jet magnetic 
field configuration. 

Jets that expand into a pre-existing, evacuated funnel are 
of different nature. They are free to accelerate to super fast- 
magnetosonic velocities, as illustrated in Fig. 1(b). The tip of 
the latter jets cannot communicate backward, and the jet mate¬ 
rial behaves as if it were part of an infinite jet. In these jets, the 
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toroidal and poloidal fields can be in equipartition in the fluid frame 
(Tchekhovskoy et al. 2008; Lyubarsky 2009, 2011). We term these 
jets headless jets to distinguish them from the headed jets discussed 
above. 

In this work we address two important questions in the physics 
of magnetised GRB jets: the stability of the jet as it propagates in 
the star and the dissipation of its magnetic energy. We present the 
results from a numerical study of the formation and propagation of 
relativistic Poynting flux dominated jets in a dense medium typical 
for a GRB progenitor star. The jets are formed by the rotation of 
a magnetised neutron star (NS) with a super-strong magnetic field, 
or a magnetar, in the same manner as described earlier. Thus, the 
magnetic field configuration is generated self-consistently, without 
any ad hoc assumptions about the ratio between B$ and B p . This 
allows us to study jet stability from first principles. Note that even 
though we focus on GRB jets, our results are general and apply to 
relativistic jets in other astrophysical systems, such as AGN jets. 

We start by describing our numerical scheme and the problem 
setup (Section 2) and reviewing the main aspects of kink instabil¬ 
ity in magnetically-dominated jets (Section 3). We then discuss the 
difference between headless jets, which are launched into a pre¬ 
existing funnel (Section 4), and headed jets, which drill their way 
though a dense ambient medium (Section 5). We show that while 
headless jets are relatively stable to kink modes, headed jets are 
kink-unstable. In Section 6, we analyse the results of our 3D jet 
simulations and discuss the effect that the kink instability has on the 
jets and the limitations of our work. We then construct an analytic 
model that describes the jet properties (Section 7). In Section 8, we 
discuss the astrophysical implications for GRB and AGN jets. We 
finish by comparing our results to previous works done on the sub¬ 
ject (Section 9) and conclude (Section 10). Throughout this paper, 
we use both spherical (r, 6 ,0) and cylindrical (. R , 0, z) coordinates. 


2 NUMERICAL SCHEME, UNITS, AND PROBLEM 
SETUP 

We carry out our simulations using the HARM code (Gammie et al 
2003), with recent improvements (McKinney 2006; Tchekhovskoy 
et al. 2007, 2009; McKinney & Blandford 2009, Tchekhovskoy et 
al. 2011). This is a static mesh, 3D, general relativistic magnetohy¬ 
drodynamic (GRMHD) code capable of following the evolution of 
high magnetisation flows while conserving mass, energy and mo¬ 
mentum to machine precision. An important advantage of the code 
is its ability to use curved grids, as we explain below. 

We run both 2D (axisymmetric) and 3D simulations on a grid 
that spans a range (r in , r out ) X (0,/r) X (0, 2n) in spherical polar co¬ 
ordinates, (r, 6 ,0). We set the boundary conditions to be reflecting 
at the poles (6 = 0,7r), free streaming at the outer radial boundary 
(r = r out ), and periodic in the 0-direction. The radial grid is uni¬ 
formly spaced in log r out to the radius r br , beyond which the grid 
becomes progressively sparse. Tables 1 and 2 summarise the pa¬ 
rameters of different models we use in this work, and we give more 
details in the corresponding sections. 

At the inner radial boundary (r = r in ), we place a perfectly 
conducting sphere carrying a monopole magnetic field. At the be¬ 
ginning of the simulation, the sphere is instantaneously spun up to 
a constant angular frequency. This generates Poynting flux domi¬ 
nated outflow that emanates from the surface of the sphere. We spin 
the inner sphere at a rather large angular frequency Q ~ 0.5c/r in so 
that the light cylinder, the cylindrical radius at which co-rotation 


velocity equals the speed of light, 


is a few times r in . This allows us to speed up the simulation while 
keeping the physics the same (Komissarov et al. 2007). Outside 
the sphere we set a static background density that resembles the 
interior of a star. The ambient density is a continuous single power- 
law distribution p a = r~ a that extends between an inner “bubble” of 
radius r b ^ r in and the outer radius of the grid. We will focus in this 
work on a - 2.5, a value of power-law index that is characteristic of 
Wolf-Rayet stars (e.g. Heger et al. 2000). We assume cold medium, 
ignore gravity, 2 and do not impart any perturbations to the initial 
conditions. 

Numerical MHD codes cannot evolve vacuum of density. 
Therefore, we set the boundary condition at r - r in and the initial 
density p in the region between r in and r b by demanding that the ter¬ 
minal Lorentz factor of the flow, which is controlled by the initial 
magnetisation of the flow, is = cr 0 = b 1 1Anpc 1 = 25-50. This 
value is high enough (» 1) to lead to a relativistically-magnetised 
outflow and at the same time low enough to keep the numerical 
noise to a minimum. 

In the context of GRBs, we normalise all length scales by R L = 
10 7 cm, the light cylinder radius of a 2 ms NS; the magnetic field 
strength by B L = 10 13 G, the magnetic field on the light cylinder; 
and all energy and density terms by B 2 L /4n. Time is measured in 
the units of light crossing time of the light cylinder, R h /c = 0.33 
ms. In these units, the dimensionless 3-velocity ft = v/c = z/t. The 
fiducial value of density, p a LM = 4500, corresponds to jets of power 
7 x 10 49 erg s -1 propagating in a 10M© star. In the context of AGNs, 
one can use a light cylinder of 1.5 x 10 13 cm and a magnetic field of 
B l ^ 10 5 G that correspond to a maximally rotating supermassive 
black hole (BH) of mass 10 8 M Q . The fiducial density in these units 
corresponds to ~ 10 -7 g cm -3 at the light cylinder (see Table 2). 
Unless mentioned otherwise, we will give physical units in the text 
in the context of GRB jets. 

Table 1 summarises the parameters of the numerical grid we 
used to carry out the simulations, while Table 2 lists the physi¬ 
cal parameters of the setup, giving them in dimensionless units as 
well as physical units specialised to GRBs and AGNs. The essen¬ 
tial parameters that control the jet stability are given in the first 
two columns of Table 2: (i) the dimensionless magnetic field en¬ 
ergy density (in units of ambient density), oc 5 2 /p L c 2 , which is a 
proxy for jet power Lj, and (ii) the power-law index of the ambient 
density distribution, a. The third parameter, oc Lj / B[Rl, reflects the 
fraction of the stellar surface area that launches the outflow, which 
is slightly different in 2D and 3D simulations, as we discuss be¬ 
low. The fourth parameter, cr 0 , sets the initial magnetisation of the 
outflow and does not influence the simulation outcome so long as 
(T 0 » 1: it affects the results only at larger distances than those 
simulated here (e.g., outside the star). 

2.1 2D simulation setup 

We carry out our 2D simulations on a spherical polar grid that is 
modified to concentrate grid cells toward the polar axis, where the 

2 Whereas gravity and ambient pressure may be important close to the jet 
launching point, they become less important at larger distances that are the 
main focus of this work. For instance, the free-fall time scale of the outer 
parts of the star is of order of 100 s, which is much longer than the time it 
takes for the jet to pierce the star, as we will show. 
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Figure 2. The shape of the grid lines used in our 2D (axisymmetric) simu¬ 
lations (for clarity, we show every fourth grid line). The grid lines collimate 
toward the rotational axis (the z-axis) thereby concentrating the resolution 
in the body of the jet. In the 3D simulations, the collimation is performed 
toward the x-axis to avoid the interaction of the jets with the coordinate 
singularity (see Sec. 2.2), with the 6 and 0 angles of the grid lines deformed 
independently, in the same way 0-lines are deformed in the Figure. 


jets form, by deforming the radial grid lines into parabolae, as seen 
in Fig. 2. We direct the rotational axis along the polar axis (6 = 0) 
and use a resolution of 1024 x 256 x 1 in the r-, 6- and 0- direc¬ 
tions, respectively (see Table 1). At this resolution the light cylinder 
(R l = 2r in ) is resolved by 5 cells, at an altitude |z| = 1000r in , and 
the jet width by 40 cells. We ran the simulations with both dipole 
and monopole magnetic field geometries and verified that in both 
cases the resultant jet properties are qualitatively the same. This 
is because in both cases the jet field lines stretch from the central 
source to the jet head and return back to the source (see Section 
5). The only essential difference is that in the dipole case there is a 
region of closed field lines around the equator of the central object, 
which modify the conditions near the base of the jet (Tchekhovskoy 
et al. 2015). We verified that for an aligned rotator, this effect does 
not substantially affect the properties of the jet at higher altitudes. 
We will therefore restrict ourselves to the monopole field. This 
choice corresponds to the assumption of a large-scale magnetic flux 
threading the progenitor star (Bucciantini et al. 2009). 


Table 1. Simulation grid parameters for the various models we present in 
this paper. See Sec. 2 for the description of the parameters in the table. 


Model 

name 

Resolution 

(Nr xN e x Nf) 

rin 

Rl 

r b 

Rl 

Rl 

r out 

~Rl 

M2 

1024 x 256 x 1 

0.5 

2 

500 

5000 

M2Cyl 

1024 x 256 x 1 

0.5 

2 

500 

5000 

M3 

128 x 96 x 192 

0.8 

0.8 

800 

8000 

M3HR 

256 x 192 x 384 

0.8 

0.8 

800 

160 

M3LP 

128 x 96 x 192 

0.8 

0.8 

800 

8000 


* The radius beyond which the radial grid becomes 
progressively sparse. 


2.2 3D simulation setup 

In the 3D runs we use a slightly modified configuration of grid 
lines and magnetic field topology. The existence of a polar coor¬ 
dinate singularity in the grid along the z-axis leads to numerical 
difficulties when matter and magnetic fields pass through the pole. 
As a result, magnetic field lines cannot cross the pole freely and 
get wrapped around the z-axis. This is not a problem in the 2D 
case, since all quantities have azimuthal symmetry, and no mass, 
momentum, or energy flows through the pole. However in the 3D 
case, we found that the imperfect flow through the axis leads to an 
artificial stability of the jets. To avoid this non-physical behaviour, 
we opted for reorienting the rotational axis of the central objet to 
point along the x-direction (see, e.g., Moll et al. 2008), and we col¬ 
limate the radial grid lines along the x-axis accordingly. Namely, 
the collimation is performed such that the 6- and 0-angles of the 
radial grid lines deformed independently. Note that when reporting 
the simulation results, we will use the natural coordinate system set 
by the jet orientation, with the z-axis pointing along the rotational 
axis. 

To avoid sending any fluxes through the z-axis, we spin not 
the entire star but only its northern and southern polar regions that 
are within 70° of the rotational axis. The rotation is uniform within 
50° of the axis and smoothly tapers off down to zero at 70°. We run 
3D simulations at two resolutions: a fiducial resolution of 128x96x 
192 cells in the r-, 6-, and 0-directions, respectively, and a high 
resolution, for which the resolution is doubled in all 3 dimensions 
(see Table 1). 


3 JET STABILITY AND ACCELERATION 

Magnetically-dominated flows are subject to various types of insta¬ 
bilities. In a narrow jet, the fastest growing instability is the kink 
(m = 1) instability (Begelman 1998; Lyubarskii 1999). It excites 
large scale helical motions in the jet and can lead to the dissipation 
of the magnetic energy or even the disruption of the entire jet. Here 
and below, we will use b for the magnetic field in the proper fluid 
fame and B for the magnetic field in the lab frame. 

The instability evolves on a time scale that is of the order of 
the Alfven travel time around the unstable region. Suppose a jet 
consists of an inner core, which is located at R < R 0 and dominated 
by the poloidal field, b p = b R e R +b z e z , and an outer region, which is 
dominated by the toroidal field, b^. Here, b p and b$ are the poloidal 
and toroidal field components measured in the proper frame, re¬ 
spectively. By construction, the pitch angle of the magnetic field in 
the proper frame, b p /b is of order unity at the edge of the core. 
Then, the growth timescale of the fastest growing kink mode in the 
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Table 2. The physical parameters of the various models we present in this paper. The two essential parameters describing the astrophysical system are the first 
two ones: the dimensionless magnetic field energy density (in units of ambient density) and the ambient density slope. See Sec. 2 for the description of the 
parameters in the table. Unless otherwise noted, cgs units are used. 


Model 

name 

4npLC 2 

Scale-fre 

a 

e parameters 

3ij t 

cro * 

GI 

P 

ms 

IB-s 

Rl 

caled ] 

10 49 

par. 

M 

«BH 

AGN-scaled par 

Mbh Rl 

ametei 

10 46 

rs 

PL 

B l 

B l R f 

10 7 

Mo 

10 8 Mo 

10 13 

10- 7 

M2 

4500 

2.5 

1 

50 

2 

1 

7 

10 

1 

2.5 

1.5 

2.3 

3 

M2Cyl 

10 5 

0 

1 

50 

2 

1 

7 

10 

1 

2.5 

1.5 

2.3 

3 

M3 

4500 

2.5 

0.76 

25 

2 

1 

5 

10 

1 

2.5 

1.5 

1.7 

3 

M3HR 

4500 

2.5 

0.76 

25 

2 

1 

5 

10 

1 

2.5 

1.5 

1.7 

3 

M3LP 

45000 

2.5 

0.76 

25 

2 

1 

0.5 

10 

1 

2.5 

1.5 

0.17 

3 


t The total Poynting flux from a rotating monopole field is 2B^R^c/3 
* Defined as b 2 /4nw at the base of the jet, where w = p + u + p is the enthalpy density. 
0 Lj denotes the power of a single jet. 


lab frame is (e.g. Appl et al. 2000): 

_ 2xRo7j b p 

hdnk — j * 

Va Oq 


( 2 ) 


where yj is the bulk Lorentz factor of the flow and v A is the Alfven 
velocity, 


va _ b 
c 3/b 2 + 4ttw 


(3) 


with b = (b 2 + b 2 ) l/2 being the total field strength in the fluid frame, 
w = p + u + p the plasma enthalpy, and u and p gas internal energy 
and pressure, respectively. 

From eq. (2) we see that strong toroidal fields tend to desta¬ 
bilise the jet and lead to the growth of kink instability, whereas 
strong poloidal fields tend to stabilise the jet against kinking. Typi¬ 
cally, it takes about 5-10 growth times for the instability to evolve 
to sufficiently affect the global structure of the magnetic field. 
Mizuno et al. (2012) found that the instability growth time is linked 
with the transverse profile of the fluid frame toroidal magnetic field. 
Kink modes grow faster in steep profiles e.g. b$ oc R _1 , while in 
shallower profiles the instability takes longer to grow. Relativistic 
motions increase the growth time in the lab frame even further due 
to the time dilation, as evident in eq. (2), and therefore increase the 
stability against kink modes even further. 

The fact that the growth of the kink modes requires Alfven 
waves to travel several times around the jet, implies that fluid ele¬ 
ments must be able to communicate efficiently across the jet; other¬ 
wise the instability cannot grow. In a highly magnetised flow, with 
b 2 » 4/rw, we have v A =* c (eq. 3), for which the condition for 
strong causal contact is (Lyubarskij 1992): 


rfij < i. (4) 

Equation (4) states that a fluid element with a Lorentz factor yj 
moving along a field line that makes an angle Qj with the jet axis, 
can communicate with the jet axis in a time that is shorter than the 
time it takes it to double its altitude. Only those regions in the jet 
that are strongly causally connected with the axis can become kink 
unstable (Lyubarskij 1992; Porth & Komissarov 2015). 


4 AXIALLY SYMMETRIC HEADLESS JETS: JETS 
MOVING THROUGH LOW-DENSITY MEDIUM 

It is convenient to study highly magnetised jets in the force- 
free approximation (e.g. Beskin et al. 2004; Narayan et al. 2007; 


Lyubarsky 2009; Tchekhovskoy et al. 2008; Bromberg et al. 2014). 
In this approximation, gas inertia is neglected, and the motion of 
the plasma is given by the drift velocity of the electromagnetic field 
(Thorne et al. 1986; Beskin 1997; Narayan et al. 2007), 


Ai - 


ExB 


E 


(5) 


where E is the electric field vector in the lab frame. In this sec¬ 
tion we consider headless jets that are free to move along a pre¬ 
existing evacuated funnel. In such jets, the lab-frame toroidal mag¬ 
netic field can be approximated as (Beskin 1997; Narayan et al. 
2007; Tchekhovskoy et al. 2008), 


- B (f> - L = QRB p /c, 


( 6 ) 


where Q is the angular frequency of the central object, and B p = 
(B 2 + B 2 ) 112 and B * are the lab-frame magnetic field strengths of 
the poloidal and toroidal components, respectively. Plugging eq. (6) 
into eq. (5), we obtain the Lorentz factor, 




B 2 „ + Bf 

P f 


B 2 


b l»l 

B 2 B 2 ' 


(7) 


where we approximate B p =* b p and b 2 ^ B 2 -E 2 , since the motions 
are mostly in the poloidal direction. In the limit b^ < b p we can 
approximate: 



The last approximate equality holds in the limit R » R L , in which 
yj scales linearly with the cylindrical radius, similar to the acceler¬ 
ation in a hydrodynamic jet (Piran et al. 1993). The approximation 
b(p < b p which leads to eq. (8) holds close to the jet axis, where 
the plasma maintains strong causal connection ( yfij < 1). If a jet 
is continuously collimated, the external collimating forces need to 
be communicated across the jet. This requires the jet to maintain 
a strong causal contact with the axis. Therefore, the linear growth 
of yj can only be maintained as long as yfij « (R/Rl) x 6j < 1. 
At larger values of 0j, the requirement that the external collimat¬ 
ing force should be communicated across the jet cross-section lim¬ 
its the Lorentz factor to a value given by eq. (4), i.e. yj ^ 1 /6j 
(Tchekhovskoy et al. 2008, 2009; Lyubarsky 2009). 

We now turn to the evaluation of jet stability to kink modes. 
In a steady state, there is a force balance in the transverse direc¬ 
tion, at the jet inlet, between the magnetic pressure force and the 
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toroidal and poloidal hoop stresses, expressed via the following 
force-balance equation: 


db- | bj | 

8/r dR + 4 nR + 4 nR curw 


= 0 , 


(9) 


where the first term is the magnetic pressure gradient, the second is 
the hoop stress of the toroidal field and the third is the hoop stress 
of the poloidal field, with R cmw denoting the poloidal curvature ra¬ 
dius of magnetic field lines (see also Tchekhovskoy et al. 2008). A 
jet cannot be self-collimated and needs to be pressure supported on 
the outside by some external confining force (Eichler 1993; Tomi- 
matsu 1994; Begelman & Li 1994; Beskin et al. 1998). In a jet 
that expands into a preexisting evacuated funnel, the force is dic¬ 
tated by the shape of the funnel walls. On the other hand, if the jet 
is pushing its way through an ambient medium, the confining force 
comes from a pressurised cocoon that is formed as a result of the jet 
propagation (see Fig. 1 and Sec. 5). The longitudinal (along the jet) 
pressure profile in the cocoon is relatively uniform and collimates 
the jet into a cylinder. This implies that the poloidal field lines are 
mostly along the z-direction, their curvature radius is very large, 
R cmw » R , and the third term in eq. (9) is negligible in comparison 
with the rest of the terms. 

In order to solve eq. (9), we need to know the relationship 
between b p and b Far from the jet axis, the velocity is mostly 
poloidal, and we can approximate b p - B p and 




, R <£ R h , 
,R^R l , 


( 10 ) 


where we used the eqs. (6) and (8) to connect B$ and B p . Thus, 
inside the light cylinder the magnetic field is dominated by the 
poloidal field, while outside the light cylinder the two are in 
equipartition, b p - b$. Substituting this into eq. (9), we get that 
the profile of the toroidal field follows 


b# oc 


R , R «: R l , 
R~ 1 ’ 2 ,R^>R h . 


( 11 ) 


Such a flat configuration of toroidal field outside R h , was shown 
by Mizuno et al. (2012) to suppress the growth rate of kink modes. 
Moreover, the larger Lorentz factors decrease the lab-frame growth 
rate even further and result in a largely stable jet (see eq. 2). Making 
use of eq. (10), we obtain for the poloidal field, 


bn 


R° 

R-l/2 


,RcR h , 


( 12 ) 


To illustrate the properties of headless jets, we run a 2D sim¬ 
ulation with a rotating sphere threaded with a monopole magnetic 
field. The sphere is placed in an empty cylindrical funnel with a 
radius R = 20, carved out in a very dense ambient medium (see 
Table 2), where the medium plays the role of a rigid, confining 
wall. The funnel radius is set to be R » R h to allow for rela¬ 
tivistic motions of the jet material. As we discussed in Section 1, 
the rotation of the central sphere generates toroidal field that ex¬ 
pands outward as an Alfven wave. The expansion is blocked by the 
confining medium, resulting in a buildup of toroidal tension which 
collimates the poloidal field lines toward the rotational axis (see 
Fig. lb). The field lines undergo an initial oscillatory phase and re¬ 
lax into a cylindrical shape along the evacuated funnel (see Bogo- 
valov & Tsinganos 2005 for a similar effect in a jet-sheath config¬ 
uration). Outside the light cylinder, the jet material propagates out¬ 
ward super-Alfvenically and cannot communicate backward with 
Alfven waves. Because of this, above the initial oscillatory phase, 
it behaves as part of an infinite, time-independent cylindrical jet. 


A snapshot of the cylindrical jet is shown in Figure 3, where 
we show from left to right, meridional slices of bb p , b^/bp and 
yfi , respectively (the colour schemes show log 10 of the plotted pa¬ 
rameters). Black solid lines track equal poloidal magnetic flux sur¬ 
faces and follow the poloidal field lines. The field lines, which 
stretch out from the central compact object, make up the jet. They 
turn around at the tip of the jet (which is located outside of Figure 3 
frame) and return back in a layer that surrounds the jet. The return¬ 
ing field lines carry no energy flux and they make up an inner mag¬ 
netised cocoon that separates the magnetic jet from the funnel wall. 
The jet and the inner cocoon are separated by a surface at which the 
poloidal field lines switch direction and b p has a local minimum, 
leading to a peak in the value of b^/bp. This surface marks the jet 
outer boundary and is located in Figure 3(c) at R } ^ 15. The ratio 
b^/bp is close to unity in the jet and the Lorentz factor increases 
from the jet axis toward the jet edge. Figure 4 shows the profiles of 
bf, b p and yj3 across the jet at the height z = 200. The strong b p 
with a comparable strength to b<p outside the light cylinder (R > 1) 
is clearly seen. Also evident is the flat profile of both field compo¬ 
nents, in agreement with eqn. (11), (12), and the linear growth of 
of yfi which is expected in such a case (eq. 8). The combination of 
a flat magnetic field profile with a dominant poloidal component, 
together with the high Lorentz factors renders headless jets stable 
to kink instability. 


5 AXIALLY SYMMETRIC HEADED JETS: JETS 
MOVING THROUGH A DENSE MEDIUM 

Magnetised jets propagating in a dense stellar envelope need to 
drill a hole through the star before they can emerge from it. Such 
jets are substantially different than the headless jets that propagate 
in a previously evacuated funnel and that we discussed in Sec. 4. 
When a jet propagates through a dense external medium at a speed 
that exceeds the ambient sound speed, it pushes the ambient gas 
in front of it, forming a bow shock, as illustrated in Fig. 1(a) with 
the solid blue line. The ambient gas that crosses the shock heats up 
and forms a cocoon around the jet (shown in green in Fig. 1). The 
cocoon applies pressure on the jet and collimates it, thus playing 
the role of the confining walls in the previous example. In this case, 
however the jet continuously expends energy on drilling the hole 
and pushing the ambient gas sideways. This causes the jet to slow 
down mainly at its upper part, or the jet head (shown in dark blue 
in Fig. 1), where most of the work is being made. 

The basic features of the jet, the bow shock and the cocoon 
are seen in Fig. 5 that shows a meridional slice of the density (left 
panel) and the pressure (right panel) in a 2D simulation of the jet. 
The shock is seen as a sharp jump in colour in the density (left) and 
pressure (right) panels. To simulate the jet we used a similar setup 
as in the headless jet case: a rotating sphere with a monopole field 
that is placed at the center of the grid and is surrounded by a cold 
ambient medium. However here the medium completely surrounds 
the sphere and follows a power-law density (see Table 1, model 
M2). As in the previous case the presence of the ambient medium 
causes the poloidal field lines to collimate along the rotational axis. 
The poloidal field lines, shown with black lines, extend from the 
jet base at the central object out to the jet head, at which they turn 
around and return back. Similar to the headless jet case, they form 
an inner magnetised cocoon that shields the magnetic jet from the 
non-magnetic cocoon of the shocked ambient gas. The outgoing 
magnetic flux, which makes up the jet, has an elongated, nearly 
cylindrical ‘cigar’ shape: its radius at the widest point is Rj ^ 30. 
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logio(bf) lo gio(V logjoCb^/bp) log 10 (yP) 


Figure 3. The 4 panels show, from left to right, the values of b(f,,b p , b^/bp 
and y/5 in a headless 2D jet. The unit of length is the light cylinder radius, 
Rl = 1. The jet propagates inside of a preexisting low-density cylindrical 
funnel (of radius R - 20/?l) carved out in a high-density ambient medium. 
The central object produces an equatorial outflow, which bounces off the 
funnel walls and re-collimates at z ~ 10, forming a pair of jets (we show 
only one of them). After several oscillations, the cylindrical radius of the 
jet settles to Rj ~ 13. The jet is surrounded by returning magnetic field 
lines, which fill the rest of the cavity but carry very little to no energy flux, 
as indicated in the leftmost panel by the low value of b<p there. The two 
middle panels show that in the jet the toroidal and poloidal magnetic field 
components are roughly in equipartition in the fluid frame, b ( p < b p . The 
rightmost panel shows that the jet Lorentz factor is largest at the edge of 
the jet, as expected from the analytic models (see Sec. 4). See Fig. 4 for 
the transverse profiles of quantities shown in this figure. Since the ambient 
density is extends out to infinity and the walls are rigid, the headless jet 
eventually reaches a steady state. See Figs. 5-7 to compare to headed jets 
that do not have a pre-existing funnel to propagate through and have to drill 
one themselves. 


The velocity of the jet head, /5 h - 0.5, is substantially lower than 
the velocity of the jet material which is practically c, at the time 
shown in Fig. 5, as we discuss below. 

To support the propagation of the jet head, the pressure in 
the jet increases, mainly due to the increase in the toroidal field 
strength. Another way to look at this is that the slower moving head 
blocks the free expansion of toroidal field which otherwise would 
stream upward with a super-fast magnetosonic velocity obtained in 
eq. (8). Since toroidal field keeps getting injected into the jet by the 
winding of the field lines at the bottom of the jet, the toroidal field 
accumulates in the jet, and the ratio of toroidal to poloidal field 
increases. 

The jet and the cocoon are in pressure balance, and the width 
of the jet is determined by the pressure profile of the cocoon along 
the jet-cocoon boundary. This profile is relatively flat in the longi¬ 
tudinal direction, leading to the near-cylindrical cigar-like jet shape 
seen in Fig. 5. Moreover, the cocoon expands sideways with a ve¬ 
locity that is much slower than the velocity of the jet material along 
the jet. Thus, the typical timescale for the cocoon pressure to drop 



Figure 4. The cross-sectional profiles of the drift 4-velocity (dash-dotted 
red), b p (dashed green) and b^ (solid blue) in a headless jet (model M2Cyl), 
taken at an altitude z = 200Rl ~ 2 x 10 9 cm. The fluid-frame toroidal mag¬ 
netic field increases inside the light cylinder, Rl = 1, peaks just outside of it 
and then drops off as b^ oc R -1 / 2 . The profile of b p is flat inside and slowly 
decreases outside of the light cylinder until the jet edge at Rj = 13Rl- The 
sharp decrease in b p at Rj occurs due to the flip in the direction of the field 
lines at the edge of the jet, marked by Rj. The returning field lines fill the 
rest of the cavity up to R\ c = 20Rl and form an inner cocoon that separates 
the jet from the medium outside the funnel. The poloidal and toroidal mag¬ 
netic field components are roughly in equipartition inside the jet, b p - b^. 
As a result (see eq. 8), the drift velocity y^fid scales linearly with R, and 
reaches a peak value of ydfid = 5 at R = 7Rl- The full morphology of the 
jet is shown in Fig. 3. 


by a factor of two is much longer than the timescale for a jet fluid 
element to double its height. This implies that instantaneously we 
can approximate the jet geometry as a steady-state cylindrical jet, 
and we can use the force-balance equation (9) to analyse the trans¬ 
verse profile of the magnetic field in the jet. 

Here again we need to know the relationship between and 
b p . Due to the accumulation of toroidal field in the jet, we approxi¬ 
mate » b p at R » R L , which leads to the following solution to 
eq. (9): 


be/, oc 


R , R <§c R h , 
R~ l , R » R L . 


(13) 


The excess of toroidal over poloidal field implies that the acceler¬ 
ation of the plasma cannot be very efficient. Substituting the ap¬ 
proximation » b p in eq. (7), we get that the drift 4-velocity is 
now 

7j/?j ~ (fif -^ 2 ) 1/2 ~ J > (headed jets) (14) 

where we assume here that \B^\ exceeds E by an order unity factor 
due to the accumulation of toroidal field. In this case eq. (6) which 
is still correct in the lab frame, holds approximately for the fluid- 
frame as well, namely b p ^ b^R^/R. This implies that 

, f constant , R R L , _ 

,*»* L . (15) 

Comparison with eq. (11) shows that the profiles of both compo¬ 
nents of the magnetic field in this case are steeper at R > R L than 
in the case of a headless jet. 
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Figure 5. Meridional slices through a 2D headed jet, i.e., a jet that propa¬ 
gates through a dense ambient medium (without any pre-drilled funnel), at a 
time t = 730/?l/c ~ 0.25 s in the fiducial 2D model M2. Density (left panel) 
and pressure (right panel) distributions are shown in colour on a log scale 
(see the colour bars). Solid black lines show magnetic field lines, which ex¬ 
tend from the surface of the central compact object out to the jet head and 
return back on the outside of the jet. Because the jet head moves supersoni¬ 
cally relative to the ambient sound speed, an arrowhead-shaped bow shock 
forms and appears as a sharp jump in colour in both panels. The ambient 
gas, which is heated at the shock, forms a cocoon that pressure supports 
the jet and causes it to assume a cigar-like shape. The bow shock expands 
sideways with time, and the cocoon pressure decreases. As a result the jet 
slowly widens and never reaches a true steady state (compare to headless 
jets in Fig. 3). 


To illustrate the magnetic field properties and the dynamics in 
the headed jets, we show in Fig. 6(a)-(d) the spatial distributions 
of b p, b p , bplb p , and yfi of the jet shown in Fig. 5. The colour 
bars are the same as in Fig. 3, to allow for easy comparison to the 
headless jets. Panel by panel comparison of Fig. 6 and Fig. 3 shows 
that in headed jets bp is stronger (panel a) and b p is much more 
concentrated toward the axis (panel b) in the headed jet, as expected 
from our analytic scalings of the field profiles (see eqs. 11-15). 
Therefore bplb p (panel c) is much higher in the headed jet than in 
the headless jet. The two red stripes that appear at R ^ 5 and R ^ 15 
mark regions where the poloidal field flips direction. The 4-velocity 
(panel d) is much lower in the headed jet than in the headless jet, 
as expected due to the need for the headed jet to drill through the 
ambient medium. 

Figure 7 shows the jet properties in a more quantitative way by 
plotting the ID profiles of b p , bp and the drift 4-velocity y d j3 d across 
the jet. The profiles are measured at z = 200 R L ~ 2 x 10 9 cm, 
the widest point of the headed jet where it has a similar width to 
the headless jet and a geometry that is close to cylindrical, making 
it easy to compare the two cases (see Figs. 5-6). It can be seen 
that b p peaks at R ~ R L and scales as b p oc R~ l beyond that, in 
agreement with eq. (13). This profile is steeper than for the headless 
jet, bp oc R ~ l/2 , as seen in Fig. 4. On the other hand, b p has a flat 


Figure 6. The 4 panels show, from left to right, meridional slices of bp, 
b p , bp I b p and yp of a headed jet drilling its way through a dense ambient 
medium, at a time t = 730Rl/c ~ 0.25 s in our fiducial 2D model M2. 
The colour scheme is the same as for the headless jet in Fig. 3 for ease of 
comparison. Poloidal field lines are shown as black solid lines. The edge of 
the jet, the boundary between outgoing and returning field lines, is located 
at Rj - 10-15/?l- Along the jet edge the poloidal field reaches a minimum 
and the ratio bp/b p reaches a maximum (panel c). See Fig. 7 for profiles of 
various quantities across the jet. 

core at R < R L and drops sharply outside of it to a very low value, 
b p ^ bp, unlike in the headless case for which b p ^ bp. As a result, 
the acceleration is much less efficient in headed jets, with y d j3 d < 2 
across the jet. The combination of these three properties - higher 
bp/bp ratio, steep profile of b p , and lower yj3, suggests that headed 
jets are unstable to kink modes. As we show in the next section, 
this leads to jets that are very different than the headless jets. 


6 KINK MODES IN HEADED POYNTING-DOMINATED 
JETS 

Kink instability is a 3D effect and cannot be studied with axially - 
symmetric simulations. It develops in regions that are dominated 
by the toroidal field (in the proper frame) and that maintain strong 
causal connection accross the jet (eq. 4). For this reason, collimated 
jets are ideal environments for the development of the kink insta¬ 
bility. Here we show that the instability evolves in two stages. First, 
it grows internally in the jet without affecting the overall jet mor¬ 
phology. This instability is often referred to as the internal kink. It 
converts the magnetic energy into thermal energy via magnetic re¬ 
connection. As a result, the toroidal magnetic field decays, and the 
jet finds itself in a stable configuration that inhibits further growth 
of the internal kink. However, kink modes can still grow externally 
on the periphery of the jet and perturb the entire jet body. Such an 
external kink instability grows over longer time scales and typically 
affects the outer parts of the jet, near the jet head. As we will see 
below, as a result of the external kink, the jet head wobbles and in- 
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Figure 7. The cross-sectional profiles of the drift 4-velocity (dash-dotted 
red) and the two fluid-frame magnetic field components, b p (dashed green) 
and (solid blue), in a headed 2D jet (model M2) at an altitude |z| = 
200/?l ~ 2xl0 9 cm. Here, the fluid-frame magnetic field scales as b$ oc R~ x 
and peaks around the light cylinder, = 1 • The profile of b p is flat inside 
and drops sharply at R » /?l, where the toroidal field dominates, b^ » b p . 
As a result, the drift velocity yafid remains of order unity across most of 
the jet. This is in contrast to the headless jet that accelerates to much higher 
velocities (see Fig. 4). The radii of the jet, which contains the outgoing 
magnetic field lines (Rj), the inner cocoon, which contains the returning 
field lines (R\ c ), and the outer cocoon, which is made up of shocked ambient 
gas (R c ), are marked with dotted vertical lines. The full morphology of the 
headed jets is shown in Figs. 5-6. 


creases its effective cross-section. This makes it harder for the jet 
to drill its way through the ambient medium and decreases the jet 
propagation velocity. 


6.1 The growth of internal kink in a collimated headed jet 

An outflow from the central compact object is initially over¬ 
pressured relative to the surrounding medium and expands in all 
directions. As the flow expands, its pressure drops until it becomes 
comparable to the surrounding pressure of the confining medium, 
and the outflow becomes collimated. 

Before the collimation point, the jet internal pressure exceeds 
the pressure of the cocoon, so the jet material expands freely. The 
Lorentz factor scales as yj oc R/R L (eq. 8) until yj = Or 1 . From this 
point on, the jet material loses the strong causal contact with the jet 
axis and the acceleration continues less efficiently (Tchekhovskoy 
et al. 2008). During the free expansion phase kink modes cannot 
grow, since their typical growth time is longer than the propagation 
time. Indeed, from eq. (2), it follows that in the linear acceleration 
regime 6dnk ~ 2 nR 2 /R L c > R/c ~ t. Beyond that the loss of strong 
causal contact inhibits the growth of the kink instability. 

At the collimation point, the jet pressure becomes equal to 
the pressure in the cocoon. The conditions for strong causal con¬ 
tact are remet and the kink instability can grow. If the collimation 
takes place outside of the linear acceleration regime, then at the 
collimation point we have > b p . The hoop stress, which be¬ 
comes effective and is proportional to b 2 (see eq. 9), cannot be 
counterbalanced by the poloidal pressure. This leads to a recolli- 
mation point—contraction of the jet cross section—followed by a 
bounce when the poloidal pressure becomes strong enough to re¬ 


sist the hoop stress, resulting in a nozzle-like shape (Lyubarsky 
2009). In a headless jet its cylindrical radius continues to oscillate 
above the nozzle until the magnetic field configuration relaxes to an 
equilibrium state, described by eq. (10), and seen in Fig. 3. Since 
- b p , this state is stable against internal kink modes (Lyubarsky 
2009, 2011). In a headed jet, as Fig. 6 shows, these two compo¬ 
nents are not in equilibrium with each other, and the magnetic field 
continues to be dominated by b$ above the collimation point (at 
R ^ R h ). This suggests that the jet can be unstable to the internal 
kink instability. 

To check this we carried out 3D simulations whose setup is 
similar to that of our 2D simulations, with a few exceptions, as de¬ 
scribed in Sec. 2.2. 3 For our fiducial model, M3, we use an ambient 
density that corresponds to a GRB jet with a luminosity 5 x 10 49 
erg/s propagating in a 10M o star (see Tabs. 1 and 2 for a full list of 
models). Figure 8 shows a volume rendering of the inner regions 
of a headed jet in our fiducial model M3. The left panel shows in 
colour the quantity log 10 (|V x B |), which is a proxy for the con¬ 
duction current density. The right panel shows the logarithm of 
magnetisation, defined here as cr = B 2 / Anw. The light blue lines 
trace out the magnetic field in the lab frame (B). Close to their 
footpoints, at z ~ 0, the field lines expand and have an ordered, 
toroidally-dominated configuration. The field lines are recollimated 
at z - 20 R l ^ 2 x 10 8 cm and converge to the jet axis and rebound, 
forming a nozzle-like shape. As they approach the nozzle, they be¬ 
gin to kink and lose their ordered shape. This is accompanied with 
a drop in the cr parameter and a peak in the conduction current den¬ 
sity, indicative of magnetic dissipation. The helical pattern of the 
dissipation region is characteristic of the kink instability. Above 
the dissipation region the field lines are less twisted and have a 
stronger poloidal component, which points to a dissipation that af¬ 
fects mostly the toroidal magnetic field component. 

Not all field lines lose the strong causal contact with the axis 
in the free expansion region. If the jet is collimated at an altitude 
Zcoii then strong causality is maintained across the field lines with an 
opening angle 6 < V^l/z co ii u P t0 Zcoii- If> in addition, the Lorentz 
factor along the field lines is not too high, yj < cr ^ 3 , the flow re¬ 
mains sub-fast magnetosonic. In this case it can “feel” the decel¬ 
erating material at the recollimation nozzle, above it. This has a 
similar effect as the slower moving head and leads to an increase 
in the toroidal field strength: toroidal field accumulates, the accel¬ 
eration becomes inefficient, and the region becomes kink unstable. 
The two conditions are met along field lines with an opening angle 
satisfying 


0<0 k = 



AC 

Zcoll U 


? Zcoll ^ 

, Zcoll > i?LO"o /3 ’ 


(16) 


and render the field lines kink unstable below the recollimation noz¬ 
zle. 

To study this effect in more detail and to better capture the 
small-scale structure of the dissipation region, we reran our fidu¬ 
cial 3D model M3 at twice the resolution in each dimension. We 


3 In 3D simulations, the rotational axis is directed not along the z-axis, as 
in 2D ones, but along the x-axis, to avoid the interaction of the jets with 
the coordinate singularity. Accordingly, the radial 3D grid lines collimate 
toward the x-axis. For simplicity of presentation, in the text we will still use 
the standard axis orientation, with the z-axis pointing along the rotational 
axis. 
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Figure 8. A snapshot of the central region in our fiducial 3D model M3 at 
t = 4400/?l/c ~ 1.5 s, when the jet head is at z = IOOORl ~ 10 10 cm, or 
about 10% of the stellar radius. The colour scheme in the left panel shows 
the log 10 (Vx/?), which is a tracer of conduction currents, and the right panel 
shows the log 10 (cr), which is a tracer of magnetisation. The solid lines how 
the held lines of b, the magnetic held as measured in the fluid frame and 
traced out in the lab frame. The degree of their twist rehects the degree to 
which the toroidal magnetic held component dominates in the huid frame. 
The jet is initially freely expanding sideways out to \z\ - 30/?l where it 
recollimates due to the conhning pressure of the cocoon. The free expansion 
region is characterised by a low current density (green colour) and a high 
cr (red colour; see the colour bar). Above this region, cr approaches unity 
(green colour) and the held lines become less toroidally twisted. Both of 
these rehect the dissipation of the toroidal magnetic held into heat above 
the collimation point. There the magnetic held lines converge to the axis 
and become unstable to the internal kink mode, which seen as a helical 
pattern in the jet (see Sec. 6.1). 


refer to this high-resolution simulation as model M3HR (see Ta¬ 
ble 1). Figures 9 and 10, show the longitudinal cross-sections of the 
jet along the y-z plane (left panel) and the x-z plane (right panel) 
in this high-resolution model. The jet propagates along the z-axis. 
The colour schemes are the same as in the corresponding 3D jets in 
Fig. 8. Regions where magnetic energy is dissipated are identified 
by a decrease in the cr parameter accompanied by the high conduc¬ 
tion current. Below the nozzle there is a cone of an opening angle 
0.1 < 6 < 0.17, which is kink unstable all the way down to the 
source. The jet becomes collimated at z co \\ - 17 (see Fig. 10), the 
opening angle of the unstable cone is in agreement with eq. (16) for 
our choice of the initial jet magnetisation of cr 0 = 25 (see Sec. 2). 

Most of the magnetic energy dissipation occurs at the colli¬ 
mation nozzle and above it. In this region magnetic field lines with 
different pitch angles converge and become kink unstable. The dis¬ 
sipation process that follows leads to the straightening of the field 
lines, by reconnecting out the local toroidal component and results 
in a gradual alignment of the field lines with the helical path of 



x y 


Figure 9. Meridional slices of log 10 (cr) through the x - z (left) and y-z 
(right) planes, in a 3D jet at time t = 2984Rl/c ~ 1 s in our high-resolution 
3D simulation, model M3HR. Note that the horizontal axis scale has been 
stretched by approximately an order of magnitude to clearly show the jet 
structure. The jets start out highly magnetised, with cr 0 =* 25. They recolli¬ 
mate at Zcoii - 17 Rl ~ 1.7 x 10 8 cm, converge onto the axis, and develop 
an internal kink instability, which is seen as small-scale yellow-red wig¬ 
gles in the jets. The decrease of the magnetisation (seen as the transition 
in colour from red to yellow) at the wiggles reflects magnetic energy dissi¬ 
pation via the internal kink instability. Some field lines, which have small 
enough opening angles to maintain strong causal contact across the jets, be¬ 
come unstable earlier and dissipate their energy at lower altitudes. These 
lines form a ‘wedge’ of lower cr that extends down to |z| — 10Rl ~ 10 8 cm. 


the jet axis. An example for such small angle reconnection can be 
seen in Figure 11. It shows a closeup of the field lines at the col¬ 
limation point, where most of the reconnection occurs. We trace 
out several magnetic field lines starting out with polar angles of 
6 = 0°, 20°, 30°, 50°, 70° on the rotating sphere. The left panel of 
Fig. 11 shows the field lines in an unconventional sense: it shows 
the field lines of proper frame magnetic field, b, that are traced out 
in the lab frame. These field lines are useful for determining the 
regions that are prone to becoming kink-unstable: this is because 
the growth of the kink instability depends on the ratio b^/bp, which 
is seen in the figure as the degree of azimuthal winding of the field 
lines. These field lines also give an idea of the proper angle be¬ 
tween two reconnecting field lines at the reconnection point. The 
right panel of Fig. 11 shows the conventional lab-frame field lines 
B. From the comparison of the two panels, we see that the proper 
field lines appear to be more disordered than the lab frame ones, 
and this can be a factor in encouraging their rapid reconnection. 
The volume rendering colour scheme shows the logarithm of the 
magnetisation, log 10 cr. The the high-cr jet appears as a dark sil¬ 
houette against the blue low-cr background. The colour of the field 
lines represents log 10 h , where h - w/p is the specific enthalpy. Re¬ 
gions where the field lines are turning red are the regions where the 
dissipation is taking place. The green line in each of the two panels 
indicates the polar field line, which starts at the north pole of the 
magnetar and follows the centre of the jet. As expected from the 
causality condition (16) and seen in Fig. 11, the green line begins 
to kink even before the jet converges onto the jet axis. The field 
lines with larger opening angles begin to kink only after the flow 
converges onto the axis. These outer field lines are initially strongly 
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Figure 10. Same as Fig. 9 but for meridional slices of the conduction current 
J oc log 10 (|V x B |). High (red) values of J indicate regions with particularly 
strong magnetic field twist and particularly strong dissipation. 


twisted around the central, green field line. However, as they dissi¬ 
pate their magnetic energy into heat, they gradually straighten out 
and become aligned with the green line (as seen in Fig. 11). 


6.2 External kink above the dissipation region 


The dissipation of magnetic field above the collimation point leads 
to two important outcomes. First, it relieves some of the toroidal 
pressure and decreases the ratio of b^/bp in the jet. As a result, 
Ro , the radius of the poloidal field dominated core, which con¬ 
trols the minimal growth time of the kink instability (see eq. 2), 
increases significantly. This process is depicted in Fig. 12 which 
shows logjoflfcp/Zfyl) in a cross section cut along the jet: the core 
starts out as a narrow spine at the base of the jet and gradually 
widens above the collimation point, at which most of the dissipa¬ 
tion occurs. Second, the dissipation generates thermal energy at the 
expense of toroidal magnetic field energy. In the presence of the 
thermal pressure, the transverse magnetic field profile flattens out. 
As we discussed in Sec. 3, both of these changes stabilise the jet 
against the further growth of kink modes at the jet inlet. In fact, we 
see that the internal kink saturates and the dissipation stops when 
the equipartition is reached between the thermal pressure and the 
magnetic pressure. 

To demonstrate the effect of the thermal pressure on the con¬ 
figuration of toroidal field, we return to the force balance equation 
(eq. 9) and include a thermal pressure term. In a cylindrical jet, the 
equation takes the form, 


db 2 tfj> dp t h 
8 ndR 4nR dR 


(17) 


where we neglect here the poloidal hoop stress term. We assume an 
approximate equipartition between the thermal and magnetic pres¬ 
sures, i.e. Pth = b 2 /8tt. The solution to eq. (17) inside the core 
where b^ <^b p , is b p = const. The resultant profile of b$ is 


b(b oc • 


R 

R -l/2 


, R < Rq , 
, R > R 0 . 


(18) 


This profile is similar to that of the headless jet (eqs. 11 and 12); the 


Figure 11. A zoom-in on the jet recollimation nozzle that is seen in Fig. 8 at 
z ~ 50Rl - 5 x 10 8 cm (model M3 at t - 4400/?l/c ~ 1.5 s). The left panel 
shows field lines in the fluid-frame b as traced out in the lab frame. This un¬ 
conventional representation of the magnetic field lines is useful for visually 
determining the degree of toroidal dominance of the magnetic field in the 
fluid-frame (see the text for details). The right panel shows the conventional 
field lines of the lab-fame magnetic field, B. The fluid-frame magnetic field 
lines are more disordered, which might reflect their readiness to reconnect 
and dissipate their energy. We show the field lines that originate at the sur¬ 
face of the magnetar with polar angles at their foot points equal to 0°, 20°, 
30°, 50°, and 70°. The colour scheme represents log 10 h on the field lines, 
where h is the specific enthalpy. The volume colour rendering represents 
log 10 cr. The jet has a higher cr than the confining cocoon and appears as 
a dark silhouette against the blue background. The field line, which origi¬ 
nates at the north pole of the magnetar, indicates the jet axis and is shown in 
green. The gradual straightening of the magnetic field lines along the green 
line is clearly seen in both panels and reflects the dissipation of the toroidal 
magnetic field due to the internal kink instability. 


only difference here is that outside the core, R > R 0 , the high ther¬ 
mal pressure replaces b p and flattens out the transverse profile of 
the toroidal field. Figure 13 shows the magnetic and thermal pres¬ 
sure profiles across the simulated jet, above the dissipation region. 
The profiles are measured in a snapshot of a jet in the model M3, 
at a time t = 4400 R^/c ^ 1.5 s (same as in Fig. 8) and an altitude 
z = 400i?L - 4 x 10 9 cm, which is well above the region where 
most of the dissipation takes place. We average the profiles in the 
azimuthal direction to smooth out the small axial asymmetries due 
to the external kink modes. We show the azimuthally averaged (b#) 
(solid blue), {b p ) (dashed green), (b) (solid grey) and (yjsnp th ) 
(dash-dotted red). As expected from eq. (18), the magnetic field 
configuration is split into two regions: i) an inner core, dominated 
by poloidal field with a flat profile that extends to R 0 ^ 10, ii) an 
outer region, dominated by toroidal field with a profile b$ oc R~ l/2 
that covers the outer jet, a.tR 0 < R < R^ - 36 R L ^ 3.6 x 10 8 cm, and 
extending into the inner cocoon, Rj < R < R lc ~ 60 - 6 x 10 8 cm. 
Within the jet, magnetic and thermal energies are indeed in equipar¬ 
tition, which is the assumption under which we obtained the scal¬ 
ing (18). 

The combination of a flat magnetic field profile and a high 
thermal pressure stabilises the jet. Mizuno et al. (2012) found that 
kink instability is suppressed in jets with toroidal field profiles as 
shallow as oc R~ l/ 2 , while Mignone et al. (2010) reached a simi¬ 
lar conclusion when they tested the stability of hot magnetised jets. 
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Figure 12. Similar to Figures 9 and 10 but for meridional slices of 
log 10 (bp/bf/)). Dark red colour shows the poloidally-dominated regions and 
their transverse extent reflects the value of Ro, the cylindrical radius of the 
jet core with b p ^> b^. We see that above the collimation point, \z\ > 30Rl, 
the value of Ro increases with increasing distance away from the central 
compact object. This increase reflects the conversion of the toroidal mag¬ 
netic field into heat by the internal kink mode. Notice that the top jet in the 
right panel moves off the image plane due to the external kink. 

Thus, above the dissipation region the hot jet material is stable to 
the internal kink modes and is expected to behave in a similar way 
to a hot, purely hydrodynamic flow. However, we find that kink in¬ 
stability can still grow around the jet perimeter. This is known as 
the external kink mode, which produces helical motions that de¬ 
form the entire body of the jet. This is seen in Fig. 14 as the large- 
scale bends of the (blue) jets. External kink grows on a time scale 
of the order of the time it takes an Alfven wave to travel around 
the jet perimeter. Typically it takes ~ 5 - 10 growth times for the 
global kink modes to grow substantially, and this is why most of 
the deformation occurs at the top, “older” parts of the jet, mainly 
near the jet head where the kink instability has the longest time to 
evolve. The helical motions of the kink-unstable head increase the 
effective cross-section of the jet and reduce its propagation veloc¬ 
ity. If the external kink grows to a high enough amplitude, it can 
even disrupt the jet and cause it to stall. 

6.3 Dissipation rates and convergence tests 

One of the important results of this work is the efficient dissipation 
of the magnetic field at and above the recollimation nozzle. This 
reduces the strength of the toroidal magnetic field component and 
leads to equipartition between thermal and magnetic energies. A 
similar transition in the magnetic field configuration was seen in 
the simulations of twisted magnetic flux tubes that underwhent dis¬ 
sipation via the internal kink (e.g. Hood et al. 2009; Gordovskyy & 
Browning 2011; Pinto et al. 2015). Once started, the dissipation in 
our simulation occurs over a length scale of < 10/?j, which is con¬ 
sistent with the expected growth time of the kink instability across 
the entire jet (see e.g. Fig. 9). This rate is also consistent with the 
dissipation rate found in other simulations (e.g. Mizuno et al. 2012; 
Porth & Komissarov 2015). 

To test the effect of the grid resolution on the dissipation rate, 
we reran our fiducial simulation, M3, with twice the resolution in 



Figure 13. Transverse profiles of the magnetic field components and the 
thermal pressure in the jet shown in Fig. 8 (fiducial model M3), taken at 
an altitude z = 400/?l ~ 4 x 19 9 cm, t = 4400/?l/c ~ 1.5 s, and aver¬ 
aged over the azimuthal direction to smooth out non-axisymmetric effects. 
We show (bp) (solid blue), ( b p ) (dashed green), ( b ) = ( yjb p + b 1 ^) (solid 

grey) and < ^8^p t h> (dash-dotted red). Here, the poloidal and toroidal field 
components are defined with respect to the magnetar’s rotational axis, and 
the actual jet axis might deviate from it due to the external kink mode. The 
magnetic field profile consists of a core of radius Ro ~ 10/?l - 10 8 cm, 
which is filled with a nearly uniform poloidal magnetic field, and is sur¬ 
rounded by a sheath, which is dominated by a toroidal field of a transverse 
profile, b(f, oc R~ 1/2 . The toroidal field dominated region covers the outer 
parts of the jet (10/?l < R < Rj - 36/?l - 3.6 x 10 8 cm) and the inner 
cocoon CRj < R < R[ c - 60/?l - 6 x 10 8 cm). The outer cocoon, which 
contains the shocked ambient matter and is dominated by thermal pressure, 
covers the region R[ c < R < R c =* 120/?l - 1.2 x 10 9 cm. 

each direction. We refer to this high-resolution simulation as model 
M3HR (see Table 1). Figure 15 shows the electromagnetic energy 
flux (dashed lines) and the thermal energy flux (solid lines) in the 
fiducial model M3 (blue) and its high-resolution version, model 
M3HR (green). For a proper comparison, the jets are taken to have 
the same z h - Because in the high resolution simulation the jets prop¬ 
agate at a velocity that is somewhat slower (by ~ 25%) than in our 
fiducial resolution simulation, the snapshots are taken at slightly 
different times. Even though the propagation velocities are differ¬ 
ent, the dissipation rates are similar. Moreover, the jet width in both 
cases is Rj ^ 10R L - 10 8 cm and the dissipation occurs over a 
range Az < 100R L ^ 10 9 cm that agrees with the expected kink 
growth timescale at this distance. We have not tested how the dis¬ 
sipation rate depends on the Riemann solvers used in the code (see 
e.g. O’Neill et al. 2012); this will be studied elsewhere. 

The ~ 25% difference in the head velocity comes from the 
fact that the high resolution jet is somewhat less stable to the ex¬ 
ternal kink than the lower resolution jet, and the wobbling motions 
of the head have a somewhat larger radius. As we show in Sec. 7, 
the head velocity is inverse proportional to the head radius. Thus 
the wobbling motions increase the effective jet radius by ~ 25%. 
Our simulations are therefore not fully converged in the context of 
the external kink and the head velocity, and it is desirable to make 
a proper convergence test. Such a test, however, is hard to carry 
out in global simulations like ours, due to the large computational 
resources required. An alternative is to perform a convergence test 
for the amplitude of the external kink instability in a small periodic 
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Figure 14. The morphology of the jets of different powers. The left panel 
shows a low-power jet in model M3LP at t = 10995Rl/c ~ 3.7 s. The 
right panel shows a high-power jet in model M3 at t = 4084/?l/c ~ 1.36 s. 
The high-power jet is 10 times more power than the lower-power jet and 
propagates approximately 2.7 times faster, and it reaches the same distance 
as the low-power jet in a shorter time. The large-scale bends of the jets are 
charactersistic of the external kink instability. The low-power jet is thinner 
and becomes external-kink unstable at a smaller distance. As a result, kink 
modes grow to a larger amplitude and deform it more strongly than the 
high-power jet. 


box. Mizuno et al. (2009) conducted such a test and concluded that 
when the jet is resolved by more than 20 cells across, the amplitude 
of the external kink is not sensitive to the resolution. In our case the 
fiducial resolution jet contains ~ 14 cells in the transverse direction 
at Zh = 160/? l - 1.6 x 10 9 cm, where we make the comparison be¬ 
tween the jets. In the high resolution jet the number of cells is twice 
as large. Thus it is likely that we are close to the full convergence. 

To sum up, our results are rather insensitive to the numerical 
resolution, suggesting the possibility that the dissipation rate and 
the jet velocity are determined by the underlying instability and not 
the details of the numerical scheme. Whereas the quantitative un¬ 
derstanding of the dissipation rate requires even higher resolution 
studies to probe the asymptotic scaling of the instability properties 
with resolution, the saturation of the instability at equipartition be¬ 
tween thermal and magnetic energies is a robust, important result 
of this work. 


7 AN ANALYTIC MODEL FOR A HOT MAGNETIC JET 

In the previous sections, we saw that the jet material dissipates 
about half of its magnetic energy via the internal kink instability. 
The dissipation continues until equipartition is reached between 
thermal and magnetic energies. Since the jet is hot and has a rel- 



Figure 15. The thermal energy flux (solid lines) and the Poynting flux 
(dashed lines), as a function of z. We show the fiducial model, M3, at 
t = 1885/?l/c ~ 0.6 s (thick blue lines) and a high resolution model, M3HR, 
at t = 2515/?l/c ~ 0.8 s (thin green lines). Both jets extend to the same dis¬ 
tance, Zh - 150, and have the same width, though the high resolution jet 
propagates at a velocity slower by 25% due to somewhat increased wob¬ 
bling motions of the head. Even though the jets show some differences at 
their heads, the dissipation rates at their recollimation points are essentially 
the same. This indicates that the dissipation could be driven by a large-scale 
instability and not affected by the microphysics. 


atively flat transverse pressure profile, it is stable to the internal 
kink modes (see Fig. 13 and Sec. 6.2). However, as Fig. 14 shows, 
that the jet (seen in dark blue colour) develops large-scale bends, 
characteristic of external kink modes. This means that, unlike hot 
hydrodynamic jets, hot magnetised jet can be unstable to the ex¬ 
ternal kink mode. In the extreme cases, it is conceivable that the 
external kink could even lead to a complete disruption of the jet. In 
our simulations, we find that the external kink causes the head of 
the jet to wobble sideways, increasing the effective cross-section of 
the jet and decreasing its propagation speed. 

Here, we investigate the conditions under which external kink 
modes can grow in the jet to a level where they affect jet propaga¬ 
tion. Motivated by the fact that the jet is hot (see Sec. 6.2), we make 
use of its similarity to a hydrodynamic jet and employ an analytic 
model for the propagation of a hydrodynamic jet in a medium by 
Bromberg et al. (2011). We calculate the conditions in the jet and 
check if the external kink has sufficient time to grow and consider¬ 
ably deform the jet. A reader interested in the analytic results may 
skip directly to eq. (30), which gives the growth rate, and eq. (32), 
which estimates the propagation velocity of the jet. 

Figure 16 illustrates our analytic jet model. We approximate 
the jet and the cocoon as concentric cylinders of the same height, 
Zh, and of radii R } and R c , respectively. We assume that the ambi¬ 
ent medium is cold and has a power law density profile, p a oc 
Here and below, we use the indices j, h, c and a to describe quan¬ 
tities in the jet, the jet head, the cocoon and the ambient medium, 
respectively. 

The jet and the cocoon are in pressure balance at their inter¬ 
face, R - Ry We approximate the pressure in the cocoon as uni¬ 
form, thus it can be evaluated as 



( 19 ) 
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Figure 16. A schematic representation of our analytic jet model. Right 
panel: the jet (shown in grey) is composed of an inner poloidal-field- 
dominated core of radius Ro, surrounded by a toroidal-field-dominated 
sheath of radius Rj. The jet is engulfed by a cocoon (light blue) of radius 
R c . The jet head, located at z = Zh, propagates at the velocity fihc and the 
cocoon expands laterally at the velocity J3 c c. Upper left panel: A lateral total 
pressure profile in a slice through the midpoint of a 2D jet in model M2 at 
time t = 730/?l/c - 2.4 s. Lower left panel: an analytic approximation for 
the pressure profile. We assume a flat pressure profile in the jet core and the 
cocoon, and in the sheath we approximate the pressure as p oc R~ l . 


where E c is the total energy in the cocoon, V c is the cocoon volume, 
and we assume here a relativistic thermal equation of state. 

The cocoon energy is injected at the jet head (see Sec. 5). Its 
value can be estimated as 

E c = f Lj(j3j-J3 h )dt, (20) 


where Lj is the jet luminosity, J3j and A are the velocities of the jet 
material and the head, respectively. The factor (j3j - A) measures 
the reduction in the energy flux due to the motion of the head. The 
volume of the cocoon is estimated as 




1 /*"(/■ 


fi c dt 


( 21 ) 


where A = V Pdp a c 2 is the sideways expansion velocity of the 
cocoon due to its own pressure. Combining eqs. (19), (20), and 
(21), and taking t - f dz^/fac we obtain (up to an order unity 
numerical factor): 


Pc * 


L]p a 




1/4 


( 22 ) 


where p a is evaluated at the jet head. A full calculation of the pres¬ 
sure is presented in Appendix A, where we follow the method of 
Bromberg et al. (2011). 

Figure 16 shows a schematic plot of the pressure profile in the 
jet. Using the scaling given in Sec. 6.2 and eq. (18) we get: 

( Po ,r<R 0 , 

<23> 

where we take R 0 ^ Rj/2 as the radius of the poloidal field domi¬ 
nated core. The jet pressure is connected to its luminosity via 

r* r 7 

Lj = 2n J ^—jPjyffijcRdR, (24) 


where T is the adiabatic index. An equipartition between the mag¬ 
netic and thermal pressures results in T = 3/2 (See Appendix A). 
Substituting eq. (23) into eq. (24), we get the jet pressure at Rj \ 


pj(R j) ^ 


2Lj 

9 tnjPjRjc 


(25) 


We can now use eqs. (22) and (25) to solve for jjRj from the con¬ 
dition of pressure balance at the at the edge of the jet: 


R ,n “ 


pL 

f 

A# 

,p c3 yf , 


1/6 


¥ 


(26) 


where ¥ = [/r(5 - a )(3 - a)/ 6] 1/3 is an integration constant of or¬ 
der unity (see Appendix A). 

External kink modes grow in the jet on a time scale of the 
order of the Alfven crossing time around the jet perimeter (eq. 2): 


. ~ 2 * R m 

fidnk — ? 

V A 


(27) 


where we made use of the fact that toroidal and poloidal fields are 
in equipartition at the jet edge. As we discussed above, typically it 
takes ~ 5-10x/ kink for the instability to develop to a level where the 
jet is considerably deformed (Mizuno et al. 2012). The time avail¬ 
able for the instability to grow is the time it takes a fluid element to 
reach the jet head: 


4yn — 


Zh 

c(P j -A) 


Zh_ 

c/3j’ 


(28) 


where we work in the limit where A A The jet will therefore 
be considerably kinked if 


IQ^kink ~ ZO^jTjlj 
4yn Zh 


(29) 


where 0.5 < tj < 1, and we made use of the fact that the jet is highly 
magnetised, i.e., that v A ~ c. Substituting R-p/j from eq. (26) we get 
an expression for the stability of the jet to the external kink mode: 


^ _ 1077fiQ nk 

4 


^ 20 



1/6 



(30) 


The propagation velocity of the jet is directly linked with the 
cross section of the jet head. In the limit where the head velocity is 
non-relativistic (y/A ^ 1), we can write (Matzner 2003; Bromberg 
etal. 2011): 


A - 




h 

X R lPaC 3 ’ 


(31) 


where R h is the cross sectional radius of the jet head. As long as the 
helical motions of the head, due to the kink instability, are not large 
enough to increase the cross section of the jet substantially, we can 
use eq. (26) to estimate the velocity of the head: 


A 




L A 

zlPaC 3 



(32) 


Figure 17 shows A as a function of the location of the jet head 
in models M3 and M3LP, representing a GRB jet with a typical 
luminosity and with a low luminosity, respectively. We compare the 
values of A measured from the simulations, shown with thick solid 
lines, to the analytic values due to eq. (32), shown with thin dashed 
lines. Values for the high (low) luminosity jet are shown in blue 
(red). We also show the values of A from eq. (30) in dash-dotted 
lines. Both jets are marginally stable at small z h (early times), and 
become increasingly more stable with time. 
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Figure 17. Jet head velocity /3h as a function of head position Zh for high- 
power (thick blue lines, model M3) and low-power (thin green lines, model 
M3LP) jets. The high-power jet is a factor 10 more luminous and propagates 
a factor of ~ 2 faster than the low-power one. We show the simulated values 
of ySh with solid lines, the analytic approximations according to eq. (32) 
with dashed lines, and the corresponding A parameters from eq. (30) with 
dash-dotted lines. Once A > 2, the velocity profile approaches the analytic 
approximation, indicating the the head wobbling motions are small, and the 
jet propagates similarly to a hydrodynamic jet. 


In the high-luminosity jet, the wobbling motions of the head 
are relatively large at Zh ^ 300. In this regime, J3 h deviates from 
the analytic solution and has a steeper dependence on z h . As the jet 
head propagates to higher altitudes it becomes more stable and the 
wobbling amplitude decreases. The effective cross-section of the 
jet decreases with increasing Zh» thus the head velocity increases at 
a faster rate than the analytic expectation. At z > 300 the simulated 
velocity profile is close to the analytic approximation (to within 
~ 10 %) and exhibits a similar profile, J3 h oc z^ 6 , reflecting the in¬ 
creased stability of the head to kink modes. In the low-luminosity 
jet, on the other hand, the wobbling motions remain large until the 
head reaches the altitude |z| = 800 R L =* 8 x 10 9 cm. Therefore the 
velocity remains below the analytic approximation and follows a 
steeper Zh-dependence. Based on these two cases, we suggest that 
the simulated J3 h approaches the analytic solution at A > 2. 

Note that for both high- and low-luminosity jets, at Zh ^ 60 
the head velocity decreases with increasing Zh- This decrease rep¬ 
resents early times in the simulation when the jet head just starts to 
propagate outward and the magnetic field on the surface of the NS 
is very high due to the buildup of the toroidal magnetic field when 
the jet was just formed. This leads to an increase in the Pointing 
flux which pushes the head at a greater force. As the jet begins to 
propagate, the excess of B$ and the Poynting flux decrease. 

Equations (30) and (32) depend on the inverse of (p a zjj). Since 
this expression decreases in density profiles steeper than z -2 , in 
such profiles the jet accelerates and becomes progressively more 
stable to kink modes with increasing distance. In contrast, in den¬ 
sity profiles shallower than z -2 , the jet would decelerate and be¬ 
come less stable with the increasing distance from the central com¬ 
pact object. In this case, the propagation velocity would be slower 
than the one obtained in eq. (32), since the wobbling of the jet head 
would increase the effective cross-section of the jet. Therefore, the 
jet would have to push a larger amount of gas to propagate. In such 


a case, we expect that there would be some altitude at which the 
head velocity will be comparable to the expansion velocity of the 
cocoon, a/ p c / p a - When this happens, the cocoon would begin to 
overtake the jet, and we expect that the jet would stall. 


8 ASTROPHYSICAL APPLICATIONS 


We now consider our analytic model for the jet propagation in the 
context of two systems that are expected to host magnetised rela¬ 
tivistic jets: GRBs and AGN. In the GRB case, we take a jet with a 
typical power of 3 X 10 49 erg/s (e.g. Guetta et al. 2005), that prop¬ 
agates in a fiducial Wolf-Rayet star of mass M* = 10M 0 , radius 
r * = 10 11 cm, and density profile p a oc r ~ a , with a - 2.5. The 
stability parameter (eq. 30) of the jet is 


AoRBfeh) - 2.6 


Zh 


(or-2)/6 


^ 49.5 r *, 11 


1/6 


(5 - a ) 2 (3 - a) 
3.125 


1/6 




(^*,34.37o3 ) 

(33) 

where 77 ~ 0.5-1 (see eq. 29) and we take ~ 1. Here, we 

used the notation A x = A x 10*. The inefficient acceleration of the 
jet material below the head implies that yj changes only weakly as 
a function of the position of the jet head and remains approximately 
2 over many decades. Eq. (33), therefore, gives A G rb ~ 2 . 6 , which 
suggests that a GRB jet is expected to be only marginally deformed 
by the external kink instability as it propagates through the star, and 
that the jet head velocity should not deviate much from the analytic 
approximation (31) for /3 h . The weak dependence of A on Zh implies 
that even close to the light cylinder, at Zh = 10 7 cm, A G rb remains 
larger than 1 , thus the jet is expected to be marginally stable also 
in the stellar interior. However, as we see from the simulations, the 
marginal stability implies substantial wobbling motions that reduce 
its propagation velocity to below the analytic expectation (31), in 
regions with A ^ 1. By substituting the GRB parameters into eq. 
(32), we obtain an estimate of the jet head velocity, 


Ai,GRB(Zh) : 


/ K(a- 2)/3 


/ £49.5 f*. 11 >0.3 




,34.3 


1/3 


(5 - £*)(3 - af 
0.625 


-1/3 


(34) 

which as we show in Sec. 7, provides a good approximation for the 
true velocity outside r ^ 3 x 10 9 cm, or outside of the inner few 
percents of the stellar radius where A > 2. Thus, the propagation 
velocity of a typical GRB jet in the star is sub-relativistic, and well 
described by eq. (34). 

The breakout time of the jet from the star can be obtained by 
integrating dz/fac, which gives 4 ^ 10 sec for the fiducial 
parameters used here. Bromberg et al. (2012, 2014) showed that 
the distribution of GRB duration times fits a typical breakout time 
of the jet from the progenitor star of about 10 sec, in agreement 
with the breakout time we obtain here. Note that this propagation 
velocity is somewhat slower than the propagation velocity of an 
equivalent hydrodynamic jet obtained analytically (Bromberg et al. 
2011) and in 2D simulations (e.g. Morsony et al. 2007; Mizuta & 
Aloy 2009). Although recent 3D simulations of hydrodynamic jets 
indicate that the two velocities are actually similar (Lopez-Camara 
et al. 2013). 

When A < 1 the wobbling motions of the jet head due to 
the external kink are large enough to decrease J3 h significantly. As 
a result, the breakout time of the jet from the star becomes much 
longer than 10 sec. A typical core-collapse GRB engine is active for 
at least as long as the GRB lasts (Sari & Piran 1997), or a few tens 
of seconds (Kouveliotou et al. 1993). If the jet power decreases, the 
jet becomes less stable and slower, and chances are that the central 
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engine will shut off before the jet head reaches the stellar edge. In 
this case, the jet will not be able to break out of the star and will 
remain buried in it. 

Such failed jets are favourable sources for producing low- 
luminosity GRBs that could occur when the cocoon breaks out of 
the stellar surface (Kulkarni et al. 1998; MacFadyen et al. 2001; 
Tan et al. 2001; Wang et al. 2007; Waxman et al. 2007; Katz et al. 
2010; Bromberg et al. 2011; Nakar & Sari 2012; Nakar 2015). The 
minimum power below which the jet becomes kink-unstable in the 
star is obtained by setting A G rb(/*) = 1 in eq. (33), giving a (single) 
jet power of 

<35) 

The associated isotropic equivalent luminosity is ZN RB iso - 2.6 x 
10 49 (#/7°) -2 erg/sec, assuming a jet with a characteristic half¬ 
opening angle of 10°. A jet with a lower power is less likely to 
produce a regular GRB (see also Woosley & Zhang 2007). Interest¬ 
ingly, the observed distribution of GRB luminosities shows a cutoff 
at isotropic equivalent luminosities lower than ~ 3 x 10 49 erg/sec 
(e.g. Cao et al. 2011). A slow, kink-unstable jet, provides a natural 
explanation for such a cutoff. 

After the jet breaks out from the star, the jet head acceler¬ 
ates and loses causal contact with the jet base. The jet, in principle, 
should relax back to the configuration of a steady state, headless jet 
that is relatively stable to kink modes and is less likely to undergo 
internal energy dissipation, as we discussed in Sec. 4. However, at 
this point the jet is sufficiently wide so that the collimation point 
is located far outside the light cylinder. In such a case, the poloidal 
field is too weak to counterbalance the contracting force of the hoop 
stress, resulting in a large pinching of the jet at the collimation noz¬ 
zle (see Sec. 6.1). The converging field lines have a large dispersion 
of pitch angles, and they converge into a small enough region that 
internal kink modes can grow and efficiently dissipate the magnetic 
energy, even without the extra compression of toroidal field com¬ 
ing from the jet head. We have run a version of our fiducial model, 
M3, in which the jet breaks out of a “stellar surface” at r - 800R L 
and verified that the dissipation at the recollimation nozzle contin¬ 
ues long after the breakout. A full analysis of the jet at the post 
breakout stage will be done in the future. 

The jet material that emerges from the star is therefore rela- 
tivistically hot, with half of its energy in thermal pressure and half 
still locked in the form of magnetic field. This means that there is 
sufficient thermal and magnetic energy to accelerate the jet mate¬ 
rial to high Lorentz factors that are inferred from the observations 
(see, e.g., Lithwick & Sari 2001; Panaitescu & Kumar 2002). The 
jet width at that time is of the order Rj ~ 10 9 cm (eq. 26) and agrees 
to an order of magnitude with the inferred size of the emission re¬ 
gion in GRB 970828 (Pe’er et al. 2007). The site of the dissipa¬ 
tion, which takes place at about z ~ 10Rj ~ 10 10 cm, is located 
sufficiently deep in the star, so that the amount of thermal pho¬ 
tons required by observations of prompt emissions can be produced 
(Vurm et al. 2013). Note that the generation of the high energy 
non-thermal emission, seen in many GRBs, potentially requires an 
additional dissipation process to take place relatively close to the 
photosphere, which can lie outside of the star (e.g. Shemi & Piran 
1990; Meszaros & Rees 2000; Bromberg et al. 2011). This topic is 
beyond the present work. 

Our results have important implications for AGN jet classifi¬ 
cation which divides AGN radio galaxies into two types: Fanaroff- 
Riley type I and II (FRI and FRII) galaxies (Fanaroff & Riley 1974). 
FRII galaxies show stable jets that typically extend over several 


hundred kiloparsecs (e.g. Hardcastle et al. 1998; Godfrey & Sha- 
bala 2013), with a well-defined hotspot at the jet head, and a co¬ 
coon that extends from the jet head backwards. The jets in FRI 
galaxies, on the other hand, are much shorter on average, extending 
over several tens of kiloparsec (e.g. Xu et al. 2000). They emit radi¬ 
ation throughout the jet body, indicative of volumetric dissipation, 
and their cocoons usually outrun the jet heads. A similar morphol¬ 
ogy to that of the FRI galaxy jets is expected when the jet head 
becomes kink unstable and its velocity drops below the expansion 
velocity of the cocoon. As we show in Sec. 7, this can occur if the 
jet propagates in a medium with a density profile that is flatter than 
r -2 . Recent studies have shown that the density profiles in clusters 
of galaxies are relatively fiat with p a oc r~ x at the inner part of the 
cluster (r < 100 kpc) and steepen to p a oc r 2 further out (e.g. 
Newman et al. 2013). In this case, a jet with a luminosity below 


L fr 
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if erg s' 


(36) 

will have a stability parameter A < 1 at r < 100 kpc and will be¬ 
come unstable to external kink before leaving the inner fiat density 
core. The wobbling motions increase the effective jet cross-section 
of the jet and reduce the propagation velocity until it becomes com¬ 
parable to the expansion velocity of the cocoon. From this point on, 
the jet is effectively stalled. The continuous injection of energy at 
the jet head results in a prominent, massive cocoon relative to the 
jet. In addition, the increasing wobbling motions can trigger dissi¬ 
pation along the jets that leads to radiation, making the jets shine 
and appear as FRI jets. In the opposite case, when L > L FR n , the 
jets are able to break out of the core before they become unsta¬ 
ble. Once out, they accelerate in the steep density profile outside 
the galaxy core and become increasingly stable. These jets are ex¬ 
pected to have more elongated, less energetic cocoons, and are ex¬ 
pected to have well-defined working surfaces at their heads where 
most of the work is being done to push the ambient medium. Such 
a morphology is observed in FRII jets. 


9 COMPARISON WITH OTHER WORKS 

The stability of magnetised jets was studied extensively via 3D sim¬ 
ulations, both in local, periodic boxes and in global simulations. 
In many studies, the jet morphology and the magnetic field con¬ 
figuration (including the pitch angle, or toroidal to poloidal field 
strength ratio) is prescribed in an ad hoc way. In addition, the jets 
are usually assumed to be of a cylindrical shape, and some initial 
perturbations are imposed on the field lines to break the axial sym¬ 
metry and jump-start the kink instability. This suggests that in such 
studies jet stability can be sensitive to the details of the initial jet 
setup (e.g., the parameters of jet injection boundary conditions, the 
initial conditions with which the jets are initialised, etc.) making it 
difficult to draw firm conclusions about jet stability. By comparing 
the properties of our self-consistently launched jets to other works, 
we can build up intuition on what effect the particular simulation 
details have on the simulation outcome and where our results fall 
in the global landscape of simulated jet configurations. 

Mizuno et al. (2012) studied the stability of a stationary jet 
with a poloidal field dominated core surrounded by toroidal field 
dominated sheath, which extended to the edge of their grid. They 
tested several configurations of toroidal and poloidal field lines in 
the sheath which characterised by different values of b^/bp. They 
found that the instability grew faster for steep lateral field profiles 
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and higher ratios of b^/bp. When the lateral toroidal profile in the 
sheath became flatter than R~ l/ 2 , the instability grew over a longer 
timescale and eventually saturated. 

Porth & Komissarov (2015) conducted a similar study us¬ 
ing jets that contained only toroidal magnetic field with a profile 
b^ ocR~ l outside a core. They replaced the poloidal field with ther¬ 
mal pressure, which dominated the core, was sub-dominant out¬ 
side, and had a flat profile outside the core (they found that the 
inclusion of poloidal field had no noticeable stabilising effect). 
Their jets, which are moderately magnetised, became kink unsta¬ 
ble over a time scale that is consistent with the expected growth 
time, ~ lCkkink- They witnessed that the energy dissipation associ¬ 
ated with the kink instability led to an increase in the jet core size, 
similar to our findings. Since their jets were infinite, the instabil¬ 
ity could grow in the jet indefinitely, and the jets were eventually 
disrupted after ~ 40t k \ nk . They also confirmed that the jets remain 
stable if they expand sideways with a velocity that is fast enough, 
so that the conditions for strong transverse causal contact are not 
met. We find that in a time-dependent case, the sideways expansion 
of the cocoon does lead to an increase in the jet radius with time. 
However, in the environments that are relevant for GRB and AGN 
jets, the rate of the expnasion is slow, and the jet remains in strong 
causal contact. 

Studies of the energy dissipation processes in stationary, non- 
relativistic magnetised jets were conducted by Hood et al. (2009); 
Gordovskyy & Browning (2011) and Pinto et al. (2015) in the 
context of solar flares. They all considered a cylindrical jet with 
poloidal field dominated core surrounded by a toroidal field domi¬ 
nated sheath with a profile steeper than b^ oc R~ l . In their setup the 
jet is surrounded by a uniform magnetic field in the z-direction, 
which provides support for the jet from the sides, similar to the ef¬ 
fect of the cocoon in our case. In all of their simulations the internal 
kink evolved in the jet on a time scale comparable to the expected 
linear growth growth time, ~ lCh kink . The internal kink resulted in 
magnetic energy dissipation over a comparable time scale. It led to 
the dissipation of most of the toroidal component of the field, and 
to a jet with mostly a poloidal field. This result is consistent with 
the internal dissipation in our collimated jets. 

Global simulations of Poynting dominated jets were con¬ 
ducted by Mignone et al. (2013). They injected a cylindrical jet 
at the lower boundary of their grid and tracked its propagation. 
The ambient medium was composed of a spherically symmetric 
supernova gas that expands homologously into a uniform density 
medium. The jets were assumed to be cylindrical, with only toroidal 
field and thermal pressure. The thermal pressure dominated in the 
core (R 0 ~ l/3Rj), and the toroidal pressure dominated outside 
the core. They tested jets with several cr parameters and Lorentz 
factors. Their jets were deformed considerably above ~ lOf^c, as 
expected. They also found indications that at the region where mag¬ 
netic energy was dissipated, field lines became more aligned with 
the average jet velocity, indicating an effective decrease of toroidal 
field in the jet. The overall morphology of the jet and its propaga¬ 
tion velocity were consistent with our results. 

Nakamura et al. (2007) carried out global, non-relativistic 
MHD simulations of AGN jets propagating in a gravitationally 
stratified isothermal medium. They injected non-rotating, magne¬ 
tised jets with a prescribed magnetic pitch b^/bp ~ 25. They found 
evidence of both internal kink at smaller distances and external kink 
at larger distances once they perturbed the ambient medium. They 
did not find any jet heating by the internal kink. 

Guan et al. (2014) carried out relativistic MHD simulations of 
non-rotating jets injected at a distance of ~ 10 3 gravitational radii 


into a uniform background medium and followed them for about 3 
orders of magnitude in distance. The transverse radius of their in¬ 
jected jets was marginally resolved by 2-3 grid cells. They found 
that when the jet material is injected with ratios of b^/bp > 10 it 
develops internal kink modes and experiences substantial dissipa¬ 
tion of magnetic energy, similar to our work. Above the dissipation 
regions, their jets seem to be considerably more bodily kinked then 
ours. This is possibly because they propagate in a uniform density 
profile that applies a larger resistive force on the jet head than in our 
case. Similar to our work, they found that the instabilities slowed 
down the jet propagation but did not disrupt the jets. 

On the analytic front, Lyubarsky (2009, 2011) obtained steady 
state solutions to the structure of collimated Poynting dominated 
jets. He showed that when the jet is collimated far from the lin¬ 
ear acceleration regime, it is focused into a recollimation nozzle 
due to the excess force of hoop stress. Based on these solutions, 
Lyubarsky (2010, 2012) postulated that if the recollimation noz¬ 
zle is narrow enough, strong dissipation should occur which will 
effectively convert the magnetic into thermal energy. Models for 
the propagation of a Poynting dominated GRB jet in a star were 
developed by Levinson & Begelman (2013) and Bromberg et al. 
(2014). Both groups obtained similar, mildly relativistic values for 
the velocity of the jet head. However, since they ignored the wob¬ 
bling motions of the head due to the external kink instability, they 
overestimated the head velocity. Based on their findings of the jet 
width, Levinson & Begelman (2013) concluded that the magnetic 
field is dissipated at the base of the jet, leading to a hydrodynamic 
jet. Bromberg et al. (2015) reached a similar conclusion based on 
the short breakout time of the jet from the star that their model pre¬ 
dicted. Indeed, we find that efficient dissipation of magnetic energy 
takes place at the recollimation nozzle, however it affects not all, 
but about half of the magnetic energy. 


10 DISCUSSION AND CONCLUSIONS 

In this work we present the results of a global numerical 3D MHD 
stability study of relativistic Poynting flux dominated jets that prop¬ 
agate in an ambient medium. The stability of a jet largely depends 
on the way it is launched and the ambient medium it interacts with. 
In order for our simulations to have predictive power when it comes 
to jet stability, we set them up to launch the jets in the same way 
as nature does it: via the magnetised rotation of a central compact 
object. The magnetised outflow thus produced has a wide opening 
angle initially. As it interacts with the ambient medium, magnetic 
pinch forces collimate it into twin narrow oppositely directed near- 
cylindrical jets. In this way, the structure of the jets is established 
self-consistently, without any a priori assumptions about, e.g., the 
jet Lorentz factor, opening angle, or magnetic pitch angle. All of 
these quantities control jet stability (see Sec. 3) and are determined 
self-consistently by the jet launching physics. In fact, our models 
are characterised by just two basic dimensionless parameters de¬ 
scribing the astrophysical system: the energy density of the poloidal 
magnetic field of the central compact object (in units of the ambient 
gas density) and the slope of the ambient gas density distribution 
(see Table 2). 

We make a distinction between the jets that propagate in an 
empty, pre-drilled funnel and the jets that have to drill their funnel 
for themselves. The former, headless jets , propagate unobstructed, 
accelerate to high Lorentz factors, and do not easily show signs 
of instability. The latter, headed jets , have much larger pitch an¬ 
gles (reflecting a stronger toroidal magnetic field component) and 
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propagate slower due to the need to expend energy on the drilling 
through the ambient gas; both of these effects make headed jets 
unstable to current-driven, non-axisymmetric instabilities. 

As headed jets run into the ambient medium, they collimate off 
of it. In fact, they recollimate and rebound, which leads to a nozzle¬ 
like shape of the jet (see Fig. 8). We show that recollimation- 
induced-compression causes the jets to become unstable to the in¬ 
ternal kink mode. This instability appears without any perturba¬ 
tions imparted to the jets or the ambient gas and leads to efficient 
dissipation of the toroidal magnetic field component. The instabil¬ 
ity operates until the equipartition between thermal and magnetic 
energies is reached: the thermal pressure stabilises the jets against 
the further growth of the internal kink (Hood & Priest 1979), and 
the dissipation stops at plasma p = /W/Wg ~ 1. 

We find that the dissipation rate is comparable to the growth 
rate of the internal kink modes, with most of the dissipation tak¬ 
ing place close to the recollimation point. Our simulations do not 
include explicit resistivity. Encouragingly, we find that the dissipa¬ 
tion rate is rather insensitive to the numerical resolution. The rate is 
consistent with that found in resistive MHD simulations of internal 
kink instability in the context of solar flares (e.g. Hood et al. 2009; 
Gordovskyy & Browning 2011; Pinto et al. 2015) and ideal rela¬ 
tivistic MHD simulations of steady state jets in periodic boxes (e.g. 
Mizuno et al. 2012; O’Neill et al. 2012). It is therefore possible that 
in some cases the dissipation rate can be determined by the rate at 
which the macroscopic instability creates and brings together field 
lines of opposite polarity and is rather insensitive to reconnection 
microphysics. 

We stress that whereas in our simulations jet recollimation 
leads to dissipation, the nature of dissipation is different from the 
usual picture of dissipation at a recollimation shock in a hydro- 
dynamic jet (Gomez et al. 1995; Marscher et al. 2008; Bromberg 
& Levinson 2009; Nalewajko & Sikora 2009; Meier 2012; Nale- 
wajko 2012; Nalewajko & Sikora 2012; Cohen et al. 2014). Even 
highly magnetised jets can undergo a recollimation shock provided 
the jet is super-fast magnetosonic. However, not much dissipation 
occurs at the shock itself, and instead the dissipation happens due to 
the activation of the internal kink instability in the post-shock/post- 
recollimation region: at the shock, the flow velocity decreases and 
increases, both of which encourage the development of an in¬ 
ternal kink mode. In the context of AGN, such a jet recollimation 
can be responsible for bright features seen in the jets at around the 
Bondi radius: for instance, the HST-1 component in the M87 jet 
(Biretta et al. 1999; Meyer et al. 2013; Hada et al. 2015) could be 
caused by the dissipation via the internal kink mode near the rec¬ 
ollimation point where the jet starts interacting with the interstellar 
medium and the bright components emerging from and/or passing 
through HST-1 might be the dissipation sites caused by the kink 
instability. 

The jets emerge from the dissipation region with the toroidal 
magnetic field component greatly weakened and the thermal jet 
content much boosted by the dissipation. This has a stabilising ef¬ 
fect on the internal kink mode and causes the dissipation to be lo¬ 
calised near the jet recollimation point. This dissipation continues 
long after the jets break out of the star (see Sec. 6.2). Since the mag¬ 
netic energy is dissipated deep inside the star, the emerging plasma 
contains enough photons to account for the GRB emission (Vurm 
et al. 2013). However, the non-thermal emission observed in many 
GRBs likely requires additional dissipation to occur close to or be¬ 
yond the photosphere (see, e.g., Giannios 2008). This can happen if 
the plasma cools by adiabatic expansion or sub-photospheric emis¬ 
sion and regains high enough magnetisation for the internal kink 


mode to be revived. It is also conceivable that if the plasma under¬ 
goes a rapid expansion phase when exiting the star, parts of it may 
attain a low enough magnetisation (Tchekhovskoy et al. 2010b), 
and the internal shocks may become effective; however, it is un¬ 
clear if these parts carry sufficient amounts of energy to power 
the prompt emission. Additional magnetic dissipation can be ex¬ 
pected if the magnetic flux threading the jets switches polarity (e.g. 
Drenkhahn & Spruit 2002; Giannios & Spruit 2005; Narayan et al. 
2011; McKinney & Uzdensky 2012; Deng et al. 2015), either due 
to the oblique magnetic field of a central magnetar producing a 
striped wind (Metzger et al. 2011) or due to the accretion of matter 
with alternating magnetic field polarity onto the central BH (e.g., 
Tchekhovskoy & Giannios 2014; Parfrey et al. 2015). 

We find that the jet residual magnetic field content beyond the 
recollimation point is sufficient to trigger the external kink insta¬ 
bility, which leads to large-scale helical motions and bends, caus¬ 
ing the jet head to wobble (see Fig. 14). As a result, the effective 
cross-section of the jets increases and the jet propagation velocity 
decreases. Using an analytic model, we show that the stability of 
our relativistically-hot jets to the external kink mode is controlled 
by the following dimensionless parameter, 

r 


where p a is the ambient density, z h is the distance to the jet head, 
and yj is the Lorentz factor of jet material (which is larger than that 
of the jet head; see Sec. 7). When A substantially exceeds unity, 
the external kink does not have sufficient time to evolve in the jets, 
and the wobbling motions remain small. In this case, the jets prop¬ 
agate at a similar velocity to hydrodynamic jets. However, when 
A < 1, the jet heads become unstable to external kink modes. The 
jet head velocity significantly decreases and becomes much smaller 
than that of hydrodynamic jets. Using our analytic model, we show 
the stability of the jets and the propagation velocity of their heads 
are closely related to the profile of the external density. In ambient 
density profiles steeper than p oc r -2 , the jets accelerate and become 
more stable with increasing distance. In flatter profiles, however, 
the jets decelerate and become exceedingly unstable. We expect 
that eventually the head velocity becomes comparable to the ex¬ 
pansion velocity of the shocked external gas, which then outruns 
the jets. This may cause the jets to stall. 

The concept of a stalled jet can explain the difference in the 
morphology of FRI and FRII radio galaxy jets (Fanaroff & Riley 
1974). The jets of FRI galaxies usually extend to several tens of 
kpc, are engulfed in a massive cocoon, and shine along the entire 
body. FRII galaxy jets, on the other hand, are much longer on aver¬ 
age (> 100 kpc), maintain straight course and have a well-defined, 
bright, hot spot at the jet head. These differences can be explained 
via the effects of the external kink instability. In the fiat density 
profiles of galaxy cluster cores, the jets decelerate as they prop¬ 
agate outwards. Therefore, a minimal power exists above which 
they are able to break out of these cores. We estimate this minimal 
power to be Lj^ n ~ 4 x 10 45 erg s -1 , adopting typical properties 
of radio galaxies (see Sec. 8 and Tab. 2). Jets with a lower power 
become external-kink-unstable inside the core. They develop wob¬ 
bling motions which increase their effective cross-section and cause 
them to slow down until they are overrun by their cocoons that con¬ 
tinuously expand. The increased wobbling motions can also trigger 
dissipation along the jets that leads to radiation, making the jets 
shine and appear as FRI jets. A higher luminosity jet will break out 
of the flat density core before becoming highly kink-unstable and 
will continue to accelerate in the steep ambient density profile that 
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surrounds the galaxy core. Such a jet will remain stable and will 
appear as an FRII jet. 

In the GRB case, the density profile of the stellar envelope is 
so steep that the jets become exceedingly more stable as they prop¬ 
agate outward. We find that typical GRB jets are marginally stable 
and propagate at a velocity < 0.5c, similar to the velocity of hy¬ 
drodynamic jets (see Fig. 17). They break out of the star in about 
~ 10 seconds. Bromberg et al. (2015) concluded that Poynting flux 
dominated jets need to dissipate their energy close to the source in 
order for the breakout time to be consistent with observations. We 
show here that such a dissipation mechanism indeed exists near the 
recollimation point of the jet. The breakout time agrees with that in¬ 
ferred from the distribution of GRB duration times (Bromberg et al. 
2012, 2015). The similarity between the propagation velocities of 
hydrodynamic and our relativistically-hot MHD jets suggests that 
it is difficult if not impossible to distinguish between magnetic and 
hydrodynamic energy injection mechanisms in the jets based on the 
observed breakout time alone. 

We find that if the power of a (single) GRB jet drops below 
~ 10 47 erg s" 1 , which corresponds to an isotropic equivalent lumi¬ 
nosity of ~ 2.6 x 10 49 (#/7°) -2 , the jet becomes so strongly unstable 
to the external kink mode inside the star that its propagation ve¬ 
locity significantly drops (see Sec. 8). This can lead to a situation 
where the GRB engine shuts off before the jet is able to breakout 
of the star. In this case, only the cocoon breaks out of the star, and 
feasibly leads to a different type of burst, which may resemble a 
low-luminosity GRB (Kulkami et al. 1998; MacFadyen et al. 2001; 
Tan et al. 2001; Wang et al. 2007; Waxman et al. 2007; Katz et al. 
2010; Bromberg et al. 2011; Nakar & Sari 2012; Nakar 2015). 

We note that whereas our simulations are carried out for a neu¬ 
tron star as the central compact object, our results are applicable 
to the case of an accreting black hole. This is because both black 
holes and neutron stars launch the jets via the azimuthal winding 
of poloidal magnetic field lines, as we discussed previously. The 
difference is that the field line winding in the case of the neutron 
star is due to the rotation of an actual, physical surface with mag¬ 
netic field lines frozen into it, and in the black hole case it is due 
to the rotation of the space-time (Tchekhovskoy 2015). Note, how¬ 
ever, that the presence of a black hole accretion disc, which we do 
not include in this work, may change the conditions at the base of 
the jet and affect the jet structure. For instance, a disc wind may 
alter the confining pressure profile of the jet. These effects will be 
studied elsewhere. 

In this work, we only studied jet propagation for 3 orders of 
magnitude in distance. This is short by about an order of magnitude 
in distance of the actual stellar surface. Carrying out longer-term 
simulations that extend up to and beyond the stellar surface are im¬ 
portant for verifying the stability of the jet out to larger distances 
and testing our analytic models in this new regime. We have not 
studied the processes that occur in the jets upon their breakout of 
the star. The breakout can be accompanied with substantial accel¬ 
eration (Tchekhovskoy et al. 2010b) and the loss of causal contact 
across the jet (Komissarov et al. 2010), both of which can affect the 
development of the instability. The simulations of jet breakout are 
important also in order to compute a first-principles structure of the 
jet outside of the star to improve the GRB prompt and afterglow 
emission modelling. 

The main finding of this work is the efficient dissipation of 
the jet magnetic field near the recollimation point via the internal 
kink instability. Though the dissipation rate is presently difficult 
to measure quantitatively due to the inherent time variability and 
potential sensitivity to the microphysics, the outcome of the inter¬ 


nal kink instability—a hot jet with thermal and magnetic energies 
in equipartition—is robust. It was obtained by multiple groups in 
various contexts using different numerical methods. Thus, the jet 
material that emerges from the star is both relativistically-hot and 
relativistically-magnetised. We conclude that the internal kink in¬ 
stability leads to jets that have enough magnetic and thermal energy 
to accelerate to the highly relativistic Lorentz factors inferred for 
GRBs, and carry a substantial thermal energy flux, which can be 
used to power the prompt emission. 
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APPENDIX A: ANALYTIC MODEL OF THE JET 

In this Appendix we calculate the jet radius following the calcu¬ 
lation method of Bromberg et al. (2011), with some minor modi¬ 
fications to fit the specific model of the magnetised jet. As illus¬ 
trated in Fig. 16, we assume a cylinder-shaped jet surrounded by 
a cylindrical cocoon propagating in a cold medium. The jet prop¬ 
agates in the z-direction and injects energy into the cocoon at the 
outermost parts jet of the jet, or at the jet head. As a result the co¬ 
coon becomes pressurised and expands sideways into the ambient 
medium, as seen in Figure 16. The calculation here is relevant for 
an external medium density with a power-law profile: p a oc z _Q \ The 
power-law assumption implies that all the quantities in the problem 
have a power-law dependency on z as well, which makes it easier 
to integrate them. 

According to eq. (19) the cocoon pressure is 
_ E c _ f —f$h)dt 

P c 3 y r / r \ 2 ’ 

c 3 n f J3 h dt(f fi c dt) 

where J3 h is the velocity of the jet head and J3j is the velocity of the 
jet material that lies below the head and that carries the jet energy 
to the head. L } is the power, or luminosity, injected into the jet by 
the engine located at the base of the jet and the term Lj(j3^ - J3 h ) 
represents the fraction of that luminosity that reaches the jet head 
and is injected into the cocoon. /3 C is the average expansion velocity 
of the cocoon in the transverse direction. Since the cocoon expands 
into a cold gas, it drives a shock into the ambient medium. We can 
approximate the expansion velocity of the cocoon as the velocity 
of that shock: j3 c = V pdp a c 2 , where p a is the average ambient 
density along the shock’s upstream. Substituting this expression in 
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eq. (Al), we get (up to integration constants, which can be calcu¬ 
lated posteriorly): 


PaC 2 L^(fij — jffh) 


37ry3h? 2 


(A2) 


Following Matzner (2003); Bromberg et al. (2011) we define a di¬ 
mensionless parameter 


L = 


h 

nRjp a c 3 



(A3) 


where Rj is the cylindrical radius of the jet. Substituting eq. (A3) in 
eq. (A2), and defining t - Zh/fac, where z h is the distance to the jet 
head, we get (up to a constant) 


pL j 

7TC 


1/2 


-1/4,-1 


PaLjC 

2~ 


1/2 


xR] PaC 3 


1/4 


(A4) 


The cocoon maintains a pressure balance with the jet along the 
jet edge (at Rj). We can calculate Pj(Rj) the equation for the total 
jet power: 


r R i r 9 

L i = 2n J Q fT 1 p jYiPi cRdR ’ 


(A5) 


where the pj is given in eq. (23) and T is the adiabatic index. To 
calculate T we can write the total jet pressure as the sum of the 
thermal pressure and the magnetic pressure: 

Pj — (Fth — l)^th + (T mag — l)w m ag — ~ + ^mag> (A6) 

where T th (T mag ) and u th (w mag ) are the adiabatic index and the en¬ 
ergy density of the gas (magnetic field) respectively. In Sec. 6.2 we 
show that above the dissipation region the thermal and magnetic 
pressures are in equipartition, implying that u th = 3w mag . Substitut¬ 
ing that in eq. (A6) we get that 

Pj ~ 2 (^th + ^magX (A7) 

This gives an adiabatic index inside the jet of 
Pi 3 

r=—p— + 1 = ( A8 > 

Wth ' ^mag ^ 

Integrating eq. (A5) with T from (A8) we get that the jet pressure 
at R } is 


pm - 


2Lj 

Rjc' 


(A9) 


The balance between the jet pressure and the cocoon pressure 
at the jet edge implies, p c (Rj) = Pj(Rj )• Using eqn. (A4) and (A9) 
we get: 


R m = 


L >< } 

u 

(XI 

(3n (5-a)(3-a)\ 


9xPi' 

U 9 ) 


(A 10) 


where we included the numerical factor obtained from the integra¬ 
tions and defined p a as the ambient density at z h - To calculate the 
numerical factor we used the following relations (see Bromberg 
et al. 2011, appendix B): p a = f p a dV/V = ^ p a ; z h = f fat = 
^fadt and R c = f fadt = ^ fadt . 
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